Spatial variation in the biochemical and isotopic composition of corals during bleaching and recovery
CB Wall, R Ritson-Williams, BN Popp, RD Gates (Limnology and Oceanography)
January 2020
Synopsis
This R-markdown has been updated to reflect a Corrigendum https://doi.org/10.1002/lno.11662. An error in original code led to the plotting of NMDS figures with incorrect symbol color; ellipses and vectors were correct and data interpretation were unchanged.
Ocean warming threatens coral reef ecosystems by disrupting the coral-alga symbiosis and causing coral bleaching. We examined bleaching (October 2014) and 3-months post bleaching recovery (Janurary 2015) in Montipora capitata and Porites compressa that either bleached or did not bleach during a 2014 bleaching event at three reef locations (HIMB, Reef 25, Reef 44) in Kāne‘ohe Bay, O‘ahu, spanning a latitudinal (north-south) gradient of oceanic input and seawater residence.
We measured changes in total biomass (ash-free dry weight) and biomass composition (proteins, lipids, carbohydrates), photopigments (total chlorophylls (a+c2)), and nutritional modes as inferred from stable isotopes (δ13C, δ15N) of the coral host and its endosymbionts Symbiodiniaceae.
# attach map GPS data shape files--can be used in one of many plotting approaches
HI<-readOGR("data/coast_n83.shp", "coast_n83") # in meters
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
sites<-read.csv("data/environmental/Reefs_lat_long.csv", header=TRUE, na.string=NA)
# add lat, long for each site and KBay
KBay=c(21.461194, -157.81225)
HIMB=c(21.435, -157.791083)
Reef25=c(21.461194, -157.82225)
Reef44=c(21.476778, -157.833611)# Make site map figure and export
# Google now require API accounts (with billing) to run ggmap and RgoogleMaps
# see: https://github.com/dkahle/ggmap
# and: https://www.r-bloggers.com/geocoding-with-ggmap-and-the-google-api/
# one approach is to use the packages below from open source satellite images and generate a series of tile-maps, as in ArcGIS.
# see: https://github.com/MilesMcBain/slippymath
library(slippymath)
library(sf)
library(mapview)
library(purrr)
library(curl)
library(glue)
# make a lat/long boundary box
bbox =
st_bbox(
c(
xmin = -157.84,
xmax = -157.783,
ymin = 21.425,
ymax = 21.49
),
crs = st_crs("+proj=longlat +ellps=WGS84")
)
bb_tile_query(bbox) # what is the boundary box dimensions for a range of zooms
# play with zoom and max tile # using above
tile_grid = bb_to_tg(bbox, zoom = 15, max_tiles = 40)
# pmap here will pull from the open API account and compile these, saving as 'images' and outputing them to the location in your directory.
images <-
pmap(tile_grid$tiles,
function(x, y, zoom){
outfile <- glue("data/API ESRI tiles/tile_{x}_{y}.jpg")
curl_download(url = glue('https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{zoom}/{y}/{x}'),
destfile = outfile)
outfile},
zoom = tile_grid$zoom)
# you can use other websites/sources for APIs, as explained at https://github.com/MilesMcBain/slippymath
library(raster)
library(rgdal)
raster_out <- tg_composite(tile_grid, images) # combine tiles
raster_to_png(raster_out, "data/API ESRI tiles/KBay.png") # produce a tile png
# this now gives you the site map from the satellite, but lat/long info absent and mapping GPS points not yet sorted out.
#######
###########
############### USING GGMAP and USER-ACCOUNT API
# new ggmap on github works with API account, download and restart R
devtools::install_github("dkahle/ggmap", ref = "tidyup")
register_google(key = "[enter API #]") # enter user API
KBay=c(-157.81225, 21.461194) # in lon/lat, sets map center
m <- get_googlemap(center = KBay, zoom = 13 ,maptype="satellite",source="google")
ggmap(m)
# this approach works if you have set up a ggmap approach using ggplot calls.KBay2=c(21.461194, -157.81225)
map <- GetMap(center = KBay2, zoom = 13, maptype = "satellite", SCALE = 2, API_console_key="xxx")
# enter use API here at "xxx"# Make site map figure and export
#######
###########
################ Using RgoogleMaps
# KBay2=c(21.461194, -157.81225) #as lat, long
# map <- GetMap(center = KBay2, zoom = 13, maptype = "satellite", SCALE = 2, API_console_key="[enter API here]") # enter use API here
# plot and export
png(file= "figures/Fig1a_sitemap.png", height=3.5, width=3.5, units="in", res=300)
par(oma=c(3,3,0,0))
PlotOnStaticMap(map, lat=sites$latitude, lon=sites$longitude, col="white", bg="coral", pch=21, cex=1.5)
TextOnStaticMap(map, lat=sites$latitude, lon=sites$longitude, labels=c("HIMB", "Reef 25", "Reef 44"), pos=4, cex=0.6, col="white", add=TRUE)
axis(1, at = LatLon2XY.centered(map, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.4, line = 0.4, col = "white", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.6)
axis(2, at = LatLon2XY.centered(map, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.4, line = 0.4, col = "white", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.6)
par(new=T, mar=c(10,10,0,0))
plot(HI, xlim=c(-158.3, -157.61), ylim=c(21.27, 21.65), lwd=0.4, col="gray", bg="white") # Oahu
rect(-157.87, 21.39, -157.71, 21.53, lwd=0.5)
box()
dev.off() #close device and save
Figure 1a. Study site location in Kāne‘ohe Bay on the windward side of the island of O’ahu, Hawai’i
Figure 1b-c. Tagged bleached and non-bleached corals of Montipora capitata (top) and Porites compress (bottom)*
Environmental data
Environmental data was collected at three site across the study period (October 2014 - January 2015). We continuously collected temperature, light (PAR) at each reef location at 2m depth; dissolved inorganic nutrients were collected at each reef ca. twice monthly, and sedimentation was determined during one month deployments, as well as a long term dataset spanning January 2015 - January 2016.
Light and Temp
Temperature profiles from HIMB (Moke o Lo’e) were obtained from NOAA. Photosynthetically Active Radiation (PAR) was used to integrate light availability over a 24h d (i.e., Daily Light Integral, DLI). Temperature recorded at each site using loggers are calcualted as daily means, maximums, and minimums. Logger failures at Reef 25 left only data for Reef 44 and HIMB.
### hourly, long term data (June 2014 - Feb 2015) in Local Standard Time (LST)
### NOAA Moke o Lo'e station 1612480
### https://tidesandcurrents.noaa.gov/stationhome.html?id=1612480
setwd('data/environmental')
temp.long<-read.csv('KBay_temp_1612480_20140601-20150228_phys.csv')
temp.long$DATE.TIME<-as.character(temp.long$DATE.TIME)
temp.long$DATE.TIME<-as.POSIXct(temp.long$DATE.TIME, format="%m/%d/%y %H:%M")
temp.long$DATE<-as.Date(temp.long$DATE.TIME, format="%m/%d/%Y")
# temp.long[which(temp.long$DATE.TIME == "2014-10-24 00:00:00"),] # at line 3481 # date of collection
# temp.long[which(temp.long$DATE.TIME == "2015-01-14 00:00:00"),] # at line 5449 # date of collection
##### Light Loggers
## Reef 44 PAR ##
light.2014<-read.csv("light_reef44_2014.csv", header=T, na.string=NA)
light.2015<-read.csv("light_reef44_2015.csv", header=T, na.string=NA)
df<-rbind(light.2014,light.2015)
df$date<-as.Date(df$date) # corrects date format
df<-df[!(df$date < "2014-12-18"),] # start at this time
df<-df[!(df$date > "2015-02-17"),] # end at this time
R44.PAR<-df
## HIMB PAR ##
light.2014<-read.csv("light_reefHIMB_2014.csv", header=T, na.string=NA)
light.2015<-read.csv("light_reefHIMB_2015.csv", header=T, na.string=NA)
df<-rbind(light.2014,light.2015)
df$date<-as.Date(df$date) # corrects date format
df<-df[!(df$date < "2014-12-18"),] #start at this time
df<-df[!(df$date > "2015-02-17"),] # end at this time
df<-df[1:5952, ] # odd NA row at row:5953, remove this
HIMB.PAR<-df
## combine dataframes
all.PAR<-as.data.frame(c(R44.PAR[, c(1,2,4)], HIMB.PAR[4]))
colnames(all.PAR)<-c("date", "time", "PAR.R44", "PAR.HIMB")
############################
# calculate daily integrated light values for each logger
# multiply calibrated values (umol photons m-2 s-1) by 0.0864 to reach mol photons m-2 s-1
# umol/s.. /1,000,000 * 24h * 60m * 60s = 0.0864...
# since averaging over dark period here, best way is to get DLI across 24h period
df.split <- split(all.PAR, f=all.PAR$date
< as.Date("2016-10-09", format="%F")) # split df by date
df.dli<-aggregate(data.frame(R44.DLI=df.split[[1]]$PAR.R44*0.0864,
HIMB.DLI=df.split[[1]]$PAR.HIMB*0.0864),
by=list(date=df.split[[1]]$date), FUN=mean)
##### Temperature Loggers
## Reef 44 Temp ##
temp.2014<-read.csv("temp_reef44_2014.csv", header=T, na.string=NA)
temp.2015<-read.csv("temp_reef44_2015.csv", header=T, na.string=NA)
df<-rbind(temp.2014,temp.2015)
df$date<-as.Date(df$date) # corrects date format
df$date.time<-as.POSIXct(paste(df$date, df$time), format="%Y-%m-%d %H:%M")
df<-df[!(df$date < "2014-10-10"),] # start at this time
df<-df[!(df$date.time >"2014-11-19 13:00:00" & df$date.time < "2014-11-21 13:00:00"), ]
df<-df[!(df$date > "2015-02-17"),] # end at this time
R44.temp<-df
## HIMB Temp ##
temp.2014<-read.csv("temp_reefHIMB_2014.csv", header=T, na.string=NA)
temp.2015<-read.csv("temp_reefHIMB_2015.csv", header=T, na.string=NA)
df<-rbind(temp.2014,temp.2015)
df$date<-as.Date(df$date) # corrects date format
df$date.time<-as.POSIXct(paste(df$date, df$time), format="%Y-%m-%d %H:%M")
df<-df[!(df$date < "2014-10-10"),] # start at this time
df<-df[!(df$date.time >"2014-11-19 13:00:00" & df$date.time < "2014-11-21 13:00:00"), ]
df<-df[!(df$date > "2015-02-17"),] # end at this time
HIMB.temp<-df
## combine dataframes
all.Temp<-as.data.frame(c(R44.temp[, c(1,2,5,4)], HIMB.temp[4]))
all.Temp<-all.Temp[-c(12386:12389),] # removing some NAs at end of dataframe
colnames(all.Temp)<-c("date", "time", "date.time", "temp.R44", "temp.HIMB")
#### temp means ####
df.split <- split(all.Temp, f=all.Temp$date
< as.Date("2016-10-09", format="%F")) # split df by date
df.mean<-aggregate(data.frame(R44.temp=df.split[[1]]$temp.R44,
HIMB.temp=df.split[[1]]$temp.HIMB),
by=list(date=df.split[[1]]$date), FUN=mean)
df.max<-aggregate(data.frame(R44.temp=df.split[[1]]$temp.R44,
HIMB.temp=df.split[[1]]$temp.HIMB),
by=list(date=df.split[[1]]$date), FUN=max)
df.min<-aggregate(data.frame(R44.temp=df.split[[1]]$temp.R44,
HIMB.temp=df.split[[1]]$temp.HIMB),
by=list(date=df.split[[1]]$date), FUN=min)####Figures
#### light, temp all plots #### SUPPLEMENTAL FIGURE 1
setwd('figures')
par(mar=c(2,3.6,1,1.3), mgp=c(2,0.5,0))
layout(matrix(c(1,2,3,4), nrow=2, byrow=TRUE))
k=1
## long term temperature
plot(WATERTEMP~ DATE.TIME, temp.long, type="n", ylab="Temperature (°C)", ylim=c(21, 31), yaxt="n", xaxt="n", xlab="", cex.lab=0.7, cex.axis=0.7)
abline(h = 28.5, col = "gray60", lwd=1.2, lty=2)
abline(v=temp.long[2761,1], col = "goldenrod", lwd=2, lty=1)
abline(v=temp.long[4729,1], col = "goldenrod", lwd=2, lty=1)
lines(temp.long$WATERTEMP~temp.long$DATE.TIME, lwd=.3, lty=1)
legend("topleft", lty=1, col=c("goldenrod"), legend=c("collections"), lwd=1.5, bty="n", cex=0.8, x.intersp = 0.3)
par(new=T) # to plot dates, use code below
plot(WATERTEMP~ DATE, temp.long, type="n", ylab="", ylim=c(21, 31), xaxt="n", xlab="", cex.lab=0.7, cex.axis=0.7)
axis.Date(1, at=seq(min(temp.long$DATE), max(temp.long$DATE), by="2 month"), format="%b '%y", cex.lab=0.7, cex.axis=0.7)
############ PAR
plot(R44.DLI ~ date, df.dli, type="n", ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")"), sep="")), ylim=c(0, 35), xaxt="n", yaxt="n", xlab="Date", cex.lab=0.7, cex.axis=0.7)
axis.Date(1, at=seq(min(df.dli$date), max(df.dli$date), by="2 week"), format="%d-%b '%y", cex.lab=0.7, cex.axis=0.7)
axis(side=2, at=c(0, 10, 20, 30), cex.lab=0.7, cex.axis=0.7)
legend("topright", lty=1, col=c("dodger blue", "coral"), legend=c("Reef 44", "HIMB"), lwd=1.5, bty="n", x.intersp=0.4, y.intersp=1, cex=0.8, inset=c(0.06, 0.0))
with(na.omit(data.frame(date=df.dli$date, DLI=rollmean(df.dli$R44.DLI, k, fill=NA))), {
lines(date, DLI, col="dodger blue", lwd=1)
})
with(na.omit(data.frame(date=df.dli$date, DLI=rollmean(df.dli$HIMB.DLI, k, fill=NA))), {
lines(date, DLI, col="coral", lwd=1)
})
#############
# plot of calibrated temp data
# mean
plot(df.mean$date, df.mean$R44.temp, type="l", col="dodger blue", xaxt="n", xlab="", ylab="Temperature (°C)", ylim=c(21, 30), cex.lab=0.7, cex.axis=0.7, lwd=1)
axis.Date(1, at=seq(min(df.mean$date), max(df.mean$date), by="3 week"), format="%d-%b '%y", cex.lab=0.7, cex.axis=0.7)
legend("topright", lty=1, col=c("dodger blue", "coral"), legend=c("Mean", " "), lwd=1.2, bty="n", x.intersp=0.3, y.intersp=0.6, cex=0.8, inset=c(0.07, 0.01))
with(df.mean, lines(df.mean$date, df.mean$HIMB.temp, col="coral", type="l", lwd=1, xlab=NA, ylab=NA, xaxt="n", yaxt="n", ylim=c(21, 31)))
############# MAX and Min Temp
plot(df.max$date, df.max$R44.temp, type="l", col="dodger blue", xaxt="n", xlab="", ylab="Temperature (°C)", ylim=c(21, 30), cex.lab=0.7, cex.axis=0.7, cex.main=0.8, lwd=1.5)
axis.Date(1, at=seq(min(df.max$date), max(df.max$date), by="3 week"), format="%d-%b '%y", cex.lab=0.7, cex.axis=0.7)
legend("topright", lty=1, lwd=c(1.7, 1.7, 0.5, 0.5), col=c("dodger blue", "coral", "dodger blue", "coral"),
legend=c("Max", " ", "Min", " "),
bty="n", x.intersp=0.2, y.intersp=1, cex=0.8, inset=c(0.05, 0.0))
with(df.max, lines(df.max$date, df.max$HIMB.temp, col="coral", type="l", lwd=1.5, xlab=NA, ylab=NA, xaxt="n", yaxt="n", ylim=c(21, 31)))
with(df.min, lines(df.min$date, df.min$R44.temp, col="dodger blue", type="l", lty=1, lwd=0.5, xlab=NA, ylab=NA, xaxt="n", yaxt="n", ylim=c(21, 31)))
with(df.min, lines(df.min$date, df.min$HIMB.temp, col="coral", type="l", lty=1, lwd=0.5, xlab=NA, ylab=NA, xaxt="n", yaxt="n", ylim=c(21, 31)))
##### save the figure and export to directory? ####
#dev.copy(pdf, "Fig.S1.Temp.PAR.pdf", width=7, height=4)
#dev.off()Figure S1 The temperature profile from NOAA weather station at Moku o Lo’e (HIMB) (top left) with coral collection (verticle lines) and local bleaching thresholds (dashed horizontal line). DLI (top right), and daily temperature mean (bottom left), and daily temperature minimum and maximum (bottom right).
Nutrients and sedimentation
This data is from bottled samples of Kāne’ohe Bay seawater from 4 Nov 2014 - 4 Feb 2015.
Seawater (ca. 100 ml) was taken from each site and filtered through a GF/F filter and stored in acid-washed bottles and frozen at -20 °C until analyzes. Dissolved inorganic nutrients—ammonia (NH3+), nitrate + nitrite (NO3- + NO2-), phosphate (PO43-), and silicate (Si(OH)4)—were measured at University of Hawai‘i at Mānoa SOEST Laboratory for Analytical Biogeochemistry using a Seal Analytical AA3 HR nutrient autoanalyzer and expressed as μmol l-1.
Rates of sedimentation at each site were measured using sediment traps (by filtered and weighing the mass of suspended particles falling into a polyvinylchloride tube (5 cm × 42 cm) capped at the base and attached to a cinder block at each reef site at a depth of 2 m. Sedimentation rates were expressed as mg sediment d-1.
Figures
#### all plots
setwd('figures')
dates<-cbind("4-Nov '14", "18-Nov '14", "26-Nov '14", "2-Dec '14", "17-Dec '14", "23-Dec '14", "30-Dec '14", "13-Jan '15", "27-Jan '15", "4-Feb '15")
par(mar=c(2,3.6,1,1.3), mgp=c(2,0.5,0))
layout(matrix(c(1,2,3,4,5,6,6,6), nrow=2, byrow=TRUE))
# Phosphate
plot(y=Reef44$Phosphate, x=Reef44$Date, xaxt="n", type="o", xlab=NA,
ylab=expression(paste("PO"[4]^"3-" ~(mu*mol~L^-1), sep="")), ylim=c(0, 0.25), pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis.Date(1, at=seq(min(Reef44$Date), max(Reef44$Date), by ="1 mon"), format="%b '%y", cex.lab=0.8, cex.axis=0.6) # plots Reef 44
#plots Reef 25
with(Reef44, lines(Reef25$Phosphate, x=Reef25$Date, type="o", pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[2]))
#plots HIMB
with(Reef44, lines(HIMB$Phosphate, x=HIMB$Date, type="o", pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[3]))
# Ammonium
plot(y=Reef44$Ammonia, x=Reef44$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("NH"[4]^"+" ~(mu*mol~L^-1), sep="")), ylim=c(0, 2.5), pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis.Date(1, at=seq(min(Reef44$Date), max(Reef44$Date), by ="1 mon"), format="%b '%y", cex.lab=0.8, cex.axis=0.6) #plots Reef 44
#plots Reef 25
with(Reef44, lines(Reef25$Ammonia, x=Reef25$Date, type="o", pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[2]))
#plots HIMB
with(Reef44, lines(HIMB$Ammonia, x=HIMB$Date, type="o", pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[3]))
# Nitrate+Nitrite
plot(y=Reef44$N.N, x=Reef44$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("NO"[3]^"-"+"NO"[2]^"-" ~(mu*mol~L^-1), sep="")), ylim=c(0, 2), pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis.Date(1, at=seq(min(Reef44$Date), max(Reef44$Date), by ="1 mon"), format="%b '%y", cex.lab=0.8, cex.axis=0.6) #plots Reef 44
#plots Reef 25
with(Reef44, lines(Reef25$N.N, x=Reef25$Date, type="o", pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[2]))
#plots HIMB
with(Reef44, lines(HIMB$N.N, x=HIMB$Date, type="o", pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[3]))
# silicate
plot(y=Reef44$Silicate, x=Reef44$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("Si(OH)"[4],~(mu*mol~L^-1), sep="")), ylim=c(0, 25), pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis.Date(1, at=seq(min(Reef44$Date), max(Reef44$Date), by ="1 mon"), format="%b '%y", cex.lab=0.8, cex.axis=0.6) # plot Reef 44
# plots Reef 25
with(Reef44, lines(Reef25$Silicate, x=Reef25$Date, type="o", pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[2]))
# plots HIMB
with(Reef44, lines(HIMB$Silicate, x=HIMB$Date, type="o", pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[3]))
# 3 month sediment.
plot(y=sed.44$sedim..g.d, x=sed.44$Time.point, xaxt="n", type="o", xlab=NA, ylab=expression(paste("sediment"~(g~d^-1), sep="")), ylim=c(0, 0.1), pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=sed.44$Time.point, labels=dates.sed, cex.axis=0.8) #plots HIMB
#plots Reef 25
with(sed.44, lines(sed.25$sedim..g.d, x=sed.25$Time.point, type="o", pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[2]))
#plots HIMB
with(sed.44, lines(sed.HIMB$sedim..g.d, x=sed.HIMB$Time.point, type="o", pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[3]))
# year sediment
plot(y=sed.year.44$sedim..g.d, x=sed.year.44$fig.dates, xaxt="n", type="n", xlab=NA, ylab=expression(paste("sediment"~(g~d^-1), sep="")), ylim=c(0, 0.2), cex.lab=1, cex.axis=0.8)
axis.Date(1, at=seq(min(sed.year.44$fig.dates), max(sed.year.44$fig.dates), by="2 mon"), format="%b '%y", cex.lab=0.8, cex.axis=0.8)
par(new=T)
# plot Reef 44
with(sed.year.44, lines(sed.year.44$sedim..g.d, x=sed.year.44$Time.point, type="o", pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[1], cex.axis=0.8))
#plot Reef 25
with(sed.year.44, lines(sed.year.25$sedim..g.d, x=sed.year.25$Time.point, type="o", pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[2], cex.axis=0.8))
# plot HIMB
with(sed.year.44, lines(sed.year.HIMB$sedim..g.d, x=sed.year.HIMB$Time.point, type="o", pch=19, lty=2, cex=0.8, lwd=1, col=plot_colors[3], cex.axis=0.8))
legend("topleft", legend=Reefs, col=plot_colors[1:3], lwd=1, pch=19, lty=2, cex=1.1, bty="n", x.intersp=0.3, y.intersp=1.2, inset=c(0.02, 0))
##### save the figure and export to directory? ####
#dev.copy(pdf, "Fig2.DIN.sed.pdf", encod="MacRoman", height=3.5, width=7)
#dev.off()Figure 2 Dissolved inorganic nutrient concentrations and sedimentation rates from all sites through time.
Models
PO43- model and posthoc
There is a significant difference among sites, with less [PO43-] at Reef 25 vs. Reef 44.
mod.PO4<-lmer(Phosphate~Location+(1|Date), data=nutrients)
anova(mod.PO4, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Location 0.0028657 0.0014328 2 18 5.0162 0.01856 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(mod.PO4, pairwise~Location)
CLD(posthoc, Letters=letters)## Location emmean SE df lower.CL upper.CL .group
## Reef 25 0.0867 0.00764 17.7 0.0706 0.103 a
## HIMB 0.0981 0.00764 17.7 0.0820 0.114 ab
## Reef 44 0.1106 0.00764 17.7 0.0945 0.127 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05# phosphate > HIMB and Reef 44, less at Reef 25NH4+ model
No significant effects.
mod.NH4<-lmer(Ammonia~Location+(1|Date), data=nutrients)
anova(mod.NH4, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Location 0.41374 0.20687 2 18 2.0229 0.1613# no effect on ammoniumN+N model and posthoc
There is a significant differences among reefs with Reef 44 having greater [N+N].
mod.N.N<-lmer(N.N~Location+(1|Date), data=nutrients)
anova(mod.N.N, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Location 0.78463 0.39232 2 18 9.3138 0.001672 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(mod.N.N, pairwise~Location)
CLD(posthoc, Letters=letters)## Location emmean SE df lower.CL upper.CL .group
## HIMB 0.354 0.0929 17.7 0.158 0.549 a
## Reef 25 0.420 0.0929 17.7 0.225 0.616 a
## Reef 44 0.725 0.0929 17.7 0.530 0.920 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05# HIMB and Reef 25, low in N+N, Reef 44 higherH4SiO4 model
No significant effects.
mod.SiO<-lmer(Silicate~Location+(1|Date), data=nutrients)
anova(mod.SiO, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Location 10.577 5.2887 2 18 0.3286 0.7242Short-term sedimentation model
A trend for increased sedimentation through time at Reef 44 and HIMB.
mod.sed<-lmer(sedim..g.d~Location+(1|Time.point), data=sed)
anova(mod.sed, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Location 0.00084174 0.00042087 2 2 5.2212 0.1607# lsmeans(mod.sed, pairwise~Location)Long-term sedimentation model and post-hoc
Over a longer sampling interval, there is a significant effect of Site, but a post-hoc shows this may only be a weak trend, although it is in agreement with short-term data. Overall, there is a long-term effect for 2x more sedimentation at Reef 44 and HIMB compared to Reef 25.
mod.sed.year<-lmer(sedim..g.d~Location+(1|Time.point), data=sed.year)
anova(mod.sed.year, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Location 0.0062363 0.0031182 2 24 3.6666 0.04078 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-(emmeans(mod.sed.year, pairwise~Location))
CLD(posthoc, Letters=letters)## Location emmean SE df lower.CL upper.CL .group
## Reef 25 0.0335 0.00895 33.7 0.0153 0.0517 a
## HIMB 0.0597 0.00895 33.7 0.0415 0.0779 a
## Reef 44 0.0610 0.00895 33.7 0.0428 0.0792 a
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05####Models
PAR model and post-hoc A statistical difference in DLI/PAR among two reef sites (HIMB > Reef 44); data logger lost at Reef 25. Overall, HIMB PAR (DLI) on average is 4.5 mol photons m-2 d-1 greater than Reef 44.
Par.diff<-mean(df.dli$HIMB.DLI)-mean(df.dli$R44.DLI)
R44=df.dli[, 1:2]; colnames(R44)=c("Date", "PAR"); R44$Location<-"Reef 44"
HIMB=df.dli[, c(1,3)]; colnames(HIMB)=c("Date", "PAR"); HIMB$Location<-"HIMB"
PAR.mod.df<-bind_rows(R44, HIMB)
mod.PAR<-lmer(PAR~Location+(1|Date), data=PAR.mod.df)
anova(mod.PAR, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Location 636.96 636.96 1 61 130.52 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(mod.PAR, pairwise~Location)
CLD(posthoc, Letters=letters)## Location emmean SE df lower.CL upper.CL .group
## Reef 44 12.0 0.581 76.8 10.9 13.2 a
## HIMB 16.5 0.581 76.8 15.4 17.7 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05Temperature model and post-hoc
Overall there are no significant differences among temperatures (daily mean, min., max.) at the sites where loggers were maintained (Reef 44 and HIMB). Differences in temperature were small and below logger resolution/accuracy.
daily mean temperature
#### using temp means
R44=df.mean[, 1:2]; colnames(R44)=c("Date", "temp"); R44$Location<-"Reef 44"
HIMB=df.mean[, c(1,3)]; colnames(HIMB)=c("Date", "temp"); HIMB$Location<-"HIMB"
Temp.mod.mean.df<-bind_rows(R44, HIMB)
mod.mean.Temp<-lmer(temp~Location+(1|Date), data=Temp.mod.mean.df)
anova(mod.mean.Temp, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Location 0.044871 0.044871 1 129 1.7169 0.1924posthoc<-emmeans(mod.mean.Temp, pairwise~Location)
CLD(posthoc, Letters=letters)## Location emmean SE df lower.CL upper.CL .group
## Reef 44 24.9 0.143 130 24.7 25.2 a
## HIMB 25.0 0.143 130 24.7 25.3 a
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05daily maximum temperature
#### using temp maxs
R44=df.max[, 1:2]; colnames(R44)=c("Date", "temp"); R44$Location<-"Reef 44"
HIMB=df.max[, c(1,3)]; colnames(HIMB)=c("Date", "temp"); HIMB$Location<-"HIMB"
Temp.mod.max.df<-bind_rows(R44, HIMB)
mod.max.Temp<-lmer(temp~Location+(1|Date), data=Temp.mod.max.df)
anova(mod.max.Temp, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Location 0.73767 0.73767 1 129 13.134 0.0004157 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(mod.max.Temp, pairwise~Location)
CLD(posthoc, Letters=letters)## Location emmean SE df lower.CL upper.CL .group
## Reef 44 25.6 0.144 132 25.3 25.9 a
## HIMB 25.7 0.144 132 25.4 26.0 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05daily minimum temperature
#### using temp mins
R44=df.min[, 1:2]; colnames(R44)=c("Date", "temp"); R44$Location<-"Reef 44"
HIMB=df.min[, c(1,3)]; colnames(HIMB)=c("Date", "temp"); HIMB$Location<-"HIMB"
Temp.mod.min.df<-bind_rows(R44, HIMB)
mod.min.Temp<-lmer(temp~Location+(1|Date), data=Temp.mod.min.df)
anova(mod.min.Temp, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Location 0.010586 0.010586 1 129 0.163 0.6871Biological metrics
Here is the data for Three Sites (north to south: Reef 44, Reef 25, HIMB), Two Periods (October 2014 Bleaching and January 2015 Recovery), in two reef coral species (Montipora capitata and Porites compressa).
# First, response variables (*physiology*) need to be normalized, as the dataframe now shows all metrics as "per milliliter of tissue slurry". Normalize the data first, then subset the dataframe for physiological response variables. Finally, a separate dataframe will be generated for isotopic data.
setwd('data')
data<-read.csv("Biological responses_BNB.final.csv", header=T, na.string=NA)
# remove columns unnecessary for final analysis
data<-data[ , !names(data) %in% "cells.ml"]
################################################
# calculate, normalized dependent variables
################################################
# helpful shorthand
SA<-data$surface.area.cm2 # surface area in cm2
blastate<-data$blastate.ml # tissue slurry blastate in ml
###############################################
### metrics of productivity and performance ###
# ug.chlorophyll.a..cm2 == ug.chl.a.ml * blastate / SA
data$chla<- (data$ug.chl.a.ml*blastate)/SA
# ug.chlorophyll.c2..cm2 == ug.chl.c2.ml * blastate / SA
data$chlc2<- (data$ug.chl.c2.ml*blastate)/SA
# total chlorophyll (a+c2)
data$chl.tot<-data$chla+data$chlc2
###########################
#### total biomass ###
# AFDW.mg..cm2 == convert AFDW g to mg, mutiply by blastate volume, divide by cm2
data$biomass<- (data$AFDW.g.ml*1000*blastate)/SA
###########################
### biomass composition ###
# the AFDW here comes from the g..AFDW of a lyophilized aliquot of tissue used in lipid extractions
# g.protein..gdw == mg.protein.ml/1000 * blastate / SA
data$prot.gdw<- (data$mg.protein.ml/1000)/(data$AFDW.freeze.ml)
# g.lipids..gdw == in this case, this is the data directly from lipid extractions of 3 ml blastate
data$lipid.gdw<- data$g.lipids.gdw
# g.carbs..gdw == mg.carbs.ml/1000 * blastate / SA
data$carb.gdw<- (data$mg.carb.ml/1000)/(data$AFDW.freeze.ml)
# kJ.energy.content..cm2 == multiply coefficents for kJ/g macromolecule (in g/gdw) and sum together
# -23.9 kJ/g protein, -17.5 kJ/g carbs, -39.5 kJ/g lipids
kJ.gdw<-((data$prot.gdw)*23.9 + (data$carb.gdw)*17.5 + (data$lipid.gdw)*39.5)
data$energy.gdw<-kJ.gdw # energy.gdw
data<-data[!(data$Sample.ID=="127"),]
#### Subset dataframes ####
############# all the data: predictors, physiology normalized, isotope data
data.full<-data[, c(1:8, 33:41, 21:32)]
########## normalized physiological data in new dataframe ###############
phys.data<-data[, c(1:8, 33:41)]
############# isotopes predictors, physiology normalized, isotope data
iso.data<-data[, c(1:8, 33, 21:32)]
################################################
# Summarize dataframe to produce means, SE, and n values.
################################################
df<-phys.data # dataframe here is physiological data
#########################################
# make shorthand for response variables #
chl.tot<-df$chl.tot
biomass<-df$biomass
prot.gdw<-df$prot.gdw
lipid.gdw<-df$lipid.gdw
carb.gdw<-df$carb.gdw
energy.gdw<-df$energy.gdw
##########################################
# make shorthand for statistical effects #
Time.point<-df$Time.point
Period<-factor(df$Period, levels=c("Bleaching", "Recovery"))
Site<-factor(df$Site, levels=c("Reef 44", "Reef 25", "HIMB"))
Species<-df$Species
Condition<-df$Condition
df$Sample.ID<-as.factor(df$Sample.ID)
Sample.ID<-df$Sample.ID
df$Pair<-as.factor(df$Pair)
Pair<-df$Pair
Depth<-df$Depth.m
df<-df[, !names(df) %in% c("chla", "chlc2")]
######################################################
# means
chl.mean<-aggregate(chl.tot~Period + Site + Species + Condition, phys.data, mean)
AFDW.mean<-aggregate(biomass~Period + Site + Species + Condition, phys.data, mean)
prot.mean<-aggregate(prot.gdw~Period + Site + Species + Condition, phys.data, mean)
lipids.mean<-aggregate(lipid.gdw~Period + Site + Species + Condition, phys.data, mean)
carbs.mean<-aggregate(carb.gdw~Period + Site + Species + Condition, phys.data, mean)
kJ.mean<-aggregate(energy.gdw~Period + Site + Species + Condition, phys.data, mean)
# SE
chl.SE<-aggregate(chl.tot~Period + Site + Species + Condition, phys.data, std.error)
AFDW.SE<-aggregate(biomass~Period + Site + Species + Condition, phys.data, std.error)
prot.SE<-aggregate(prot.gdw~Period + Site + Species + Condition, phys.data, std.error)
lipids.SE<-aggregate(lipid.gdw~Period + Site + Species + Condition, phys.data, std.error)
carbs.SE<-aggregate(carb.gdw~Period + Site + Species + Condition, phys.data, std.error)
kJ.SE<-aggregate(energy.gdw~Period + Site + Species + Condition, phys.data, std.error)
df.summary<-cbind(chl.mean, AFDW.mean[c(0,5)], prot.mean[c(0,5)], lipids.mean[c(0,5)], carbs.mean[c(0,5)], kJ.mean[c(0,5)], chl.SE[c(0,5)], AFDW.SE[c(0,5)], prot.SE[c(0,5)], lipids.SE[c(0,5)], carbs.SE[c(0,5)], kJ.SE[c(0,5)]); colnames(df.summary)= c("Period", "Site", "Species", "Condition", "chl", "AFDW", "prot", "lipids", "carbs", "kJ", "chlSE", "AFDWSE", "protSE", "lipidsSE", "carbsSE", "kJSE")
########### sample sizes ################
# means
chl.n<-aggregate(chl.tot~Period + Site + Species + Condition, phys.data, length)
AFDW.n<-aggregate(biomass~Period + Site + Species + Condition, phys.data, length)
prot.n<-aggregate(prot.gdw~Period + Site + Species + Condition, phys.data, length)
lipids.n<-aggregate(lipid.gdw~Period + Site + Species + Condition, phys.data, length)
carbs.n<-aggregate(carb.gdw~Period + Site + Species + Condition, phys.data, length)
kJ.n<-aggregate(energy.gdw~Period + Site + Species + Condition, phys.data, length)
df.sample.size<-cbind(chl.n, AFDW.n[c(0,5)], prot.n[c(0,5)], lipids.n[c(0,5)], carbs.n[c(0,5)], kJ.n[c(0,5)])
###########
# symbionts: what does the community of M. capitata, based on Cunning et al. 2016 data
Symbs<-phys.data[(phys.data$Species=="MC"),]
# table(Symbs$Period, Symbs$Condition:Symbs$symb.comm, exclude=NULL)Multivariate analysis
NMDS
Using metaNMDS, plot NMDS plots for each species separately. For each species, produce 4 plots: 2 plots for Bleaching (Oct 2014) and 2 plots for Recovery (Jan 2015), within each period showing Site:Condition and Condition effects. Vectors (lines) show significant univariate effects.
NMDS plots for Montipora capitata
# multivariate test using adonis-- Manova, mutivariate analysis of variance
# names(data.full)
data.full[is.na(data.full$energy.gdw), ] # where are NAs? All in this column
df.manova<-data.full[!is.na(data.full$energy.gdw),] # remove rows with NA
# changed to absolute values for the MANOVA (cannot have NA or negatives)
df.manova$host.d13C<-abs(df.manova$host.d13C)
df.manova$symb.d13C<-abs(df.manova$symb.d13C)
# remove columns unnecessary for final analysis
df.manova<-df.manova[ , !names(df.manova) %in% c("prot.cm2", "AFDW.freeze", "symb.comm", "chla", "chlc2", "host.ug.N", "host.ug.C", "symb.ug.N", "symb.ug.C", "d15N.H.S", "d13C.H.S")]
Manova.data<-df.manova[, c(9:20)] # data sans factors
###################
### NMDS plot: Montipora 2014 ###
df.manova2.MC<-df.manova[(df.manova$Species=="MC" & df.manova$Period=="Bleaching"),]
NMDS<-subset(df.manova2.MC, select = c(-Time.point, -Period, -Site, -Species, -Condition, -Sample.ID, -Pair, -Depth.m))
colnames(NMDS)<-c("chl", "biomass", "proteins", "lipids", "carbs", "kJ", "d15N-host", "d13C-host", "C:N-host", "d15N-symb", "d13C-symb", "C:N-symb")
NMDS.MC<-metaMDS(NMDS, distance='bray', k=2, trymax=100)
# applies sqre root trans. Calculated Bray-Curtis distance
# iterations until solution reached (stress minimized after reconfiguration in dimesions = k)
# trymax = nuber of iterations, can help in non-convergence
# high stress? increase # of dimensions through 'k='
# stressplot(NMDS.MC)
# Shepard plot, scatter around regression between interpoints distance in final configuration
# (cont.) against their ordinal dissimilarities
# large scatter==original dissimilarites not well preserved in reduced # dimensions
factor.df<-df.manova2.MC[, c(2:8)]
# names(factor.df)
NMDS.df<-data.frame(x=NMDS.MC$point[,1], y=NMDS.MC$point[,2],
Period=as.factor(factor.df[,1]),
Site=as.factor(factor.df[,2]),
Species=as.factor(factor.df[,3]),
Condition=as.factor(factor.df[,4]),
Sample.ID=as.factor(factor.df[,5]))
NMDS.df<-NMDS.df[, c(3:7,1,2)] # reorders columns
colnames(NMDS.df)[6:7]<-c("NMDS1", "NMDS2")
# vectors for variables
vars<-envfit(NMDS.MC, NMDS, permu=999)
#scores(vars, "vectors")
##### ##### ##### ##### ##### ##### #####
##### ##### ##### ##### ##### #####
# PLOTS
### "Montipora 2014---Site:Condition"
par(mfcol=c(2,2), mar=c(4,4,1,1), pty="sq")
MDS.levels<-c("R44-B", "R44-NB", "R25-B", "R25-NB", "HIMB-B", "HIMB-NB")
MDS.colors<-c("lightpink", "firebrick2", "slategray1", "dodgerblue2","khaki", "darkorange")
NMDS.mean=aggregate(NMDS.df[,6:7],list(group=NMDS.df$Site:NMDS.df$Condition), mean)
NMDS.symbols<-c(25, 25, 23, 23, 21, 21)
groups<-NMDS.df$Site:NMDS.df$Condition
NMDS.plot<-ordiplot(NMDS.MC, type='n', main=substitute(paste(italic("Montipora capitata"))), cex.main=1, display="sites", xlim=c(-0.2, 0.2), ylim=c(-0.2, 0.2), cex.lab=0.8, cex.axis=0.8, xaxt='n', xlab="")
axis(side = 1, labels = FALSE, tck = -0.03)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
text(-0.11, 0.19, "Oct '14: Bleaching", cex=1)
points(NMDS.plot, "sites", pch=NMDS.symbols[groups], cex=1, col="black", bg=MDS.colors[groups])
with(NMDS.plot, points(NMDS.mean[,2], NMDS.mean[,3], pch=4))
ordiellipse(NMDS.plot, groups=NMDS.df$Site:NMDS.df$Condition, draw="polygon", col=MDS.colors, label=F, alpha=150, kind="se", conf=0.95)
##### ##### ##### ##### ##### ##### #####
##### ##### ##### ##### ##### #####
### "Montipora 2014---Condition"
MDS.levels<-c("Bleached", "Non-bleached")
MDS.colors<-c("lightpink", "slategray1")
NMDS.mean=aggregate(NMDS.df[,6:7],list(group=NMDS.df$Condition), mean)
NMDS.symbols<-c(21, 23)
groups2<-NMDS.df$Condition
NMDS.plot<-ordiplot(NMDS.MC, type='n', display="sites", xlim=c(-0.2, 0.2), ylim=c(-0.2, 0.2), cex.lab=0.8, cex.axis=0.8, ylab="NMDS2", xlab="NMDS1")
axis(side = 1, labels = FALSE, tck = -0.03)
axis(side = 2, labels = FALSE, tck = -0.03)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
#text(-0.25, 0.19, "b", cex=1.5)
points(NMDS.plot, "sites", pch=NMDS.symbols[groups2], cex=1, col="black", bg=MDS.colors[groups2])
with(NMDS.plot, points(NMDS.mean[,2], NMDS.mean[,3], pch=4))
ordiellipse(NMDS.plot, groups=NMDS.df$Condition, draw="polygon", col=MDS.colors, label=F, alpha=150, kind="se", conf=0.95)
par.new=T
plot(vars, p.max=0.05, cex=0.8, lwd=1)
##### ##### ##### ##### ##### ##### #####
##### ##### ##### ##### ##### #####
######### Montipora 2015 ########
df.manova2.MC<-df.manova[(df.manova$Species=="MC" & df.manova$Period=="Recovery"),]
NMDS<-subset(df.manova2.MC, select = c(-Time.point, -Period, -Site, -Species, -Condition, -Sample.ID, -Pair, -Depth.m))
colnames(NMDS)<-c("chl", "biomass", "proteins", "lipids", "carbs", "kJ", "d15N-H", "d13C-H", "C:N-H", "d15N-S", "d13C-S", "C:N-S")
NMDS.MC<-metaMDS(NMDS, distance='bray', k=2, trymax=100)
factor.df<-df.manova2.MC[, c(2:8)]; names(factor.df)
NMDS.df<-data.frame(x=NMDS.MC$point[,1], y=NMDS.MC$point[,2],
Period=as.factor(factor.df[,1]),
Site=as.factor(factor.df[,2]),
Species=as.factor(factor.df[,3]),
Condition=as.factor(factor.df[,4]),
Sample.ID=as.factor(factor.df[,5]))
NMDS.df<-NMDS.df[, c(3:7,1,2)] # reorders columns
colnames(NMDS.df)[6:7]<-c("NMDS1", "NMDS2")
# vectors for variables
vars<-envfit(NMDS.MC, NMDS, permu=999)
# scores(vars, "vectors")
##### ##### ##### ##### ##### ##### #####
##### ##### ##### ##### ##### #####
# PLOTS
### "Montipora 2015---Site:Condition"
MDS.levels<-c("R44-B", "R44-NB", "R25-B", "R25-NB", "HIMB-B", "HIMB-NB")
MDS.colors<-c("lightpink", "firebrick2", "slategray1", "dodgerblue2","khaki", "darkorange")
NMDS.mean=aggregate(NMDS.df[,6:7],list(group=NMDS.df$Site:NMDS.df$Condition), mean)
NMDS.symbols<-c(25, 25, 23, 23, 21, 21)
NMDS.plot<-ordiplot(NMDS.MC, type='n', display="sites", xlim=c(-0.2, 0.2), ylim=c(-0.2, 0.2), cex.lab=0.8, cex.axis=0.8, xaxt='n', yaxt='n', ylab="", xlab="")
axis(side = 1, labels = FALSE, tck = -0.03)
axis(side = 2, labels = FALSE, tck = -0.03)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
text(-0.11, 0.19, "Jan '15: Recovery", cex=1)
points(NMDS.plot, "sites", pch=NMDS.symbols[groups], cex=1, col="black", bg=MDS.colors[groups])
with(NMDS.plot, points(NMDS.mean[,2], NMDS.mean[,3], pch=4))
ordiellipse(NMDS.plot, groups=NMDS.df$Site:NMDS.df$Condition, draw="polygon", col=MDS.colors, label=F, alpha=150, kind="se", conf=0.95)
legend("topright", legend=MDS.levels, pch=NMDS.symbols, col="black", cex=0.85, inset=c(0.06,0), pt.bg=MDS.colors, pt.cex=1, bty="n", x.intersp = 1.2, y.intersp = 1)
##### ##### ##### ##### ##### ##### #####
##### ##### ##### ##### ##### #####
### "Montipora 2015---Condition"
### Condition
MDS.levels<-c("Bleached", "Non-bleached")
MDS.colors<-c("lightpink", "slategray1")
NMDS.mean=aggregate(NMDS.df[,6:7],list(group=NMDS.df$Condition), mean)
NMDS.symbols<-c(21, 23)
NMDS.plot<-ordiplot(NMDS.MC, type='n', display="sites", xlim=c(-0.2, 0.2), ylim=c(-0.2, 0.2), cex.lab=0.8, cex.axis=0.8, ylab="", yaxt='n')
axis(side = 2, labels = FALSE, tck = -0.03)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
#text(-0.25, 0.19, "d", cex=1.5)
points(NMDS.plot, "sites", pch=NMDS.symbols[groups2], cex=1, col="black", bg=MDS.colors[groups2])
with(NMDS.plot, points(NMDS.mean[,2], NMDS.mean[,3], pch=4))
ordiellipse(NMDS.plot, groups=NMDS.df$Condition, draw="polygon", col=MDS.colors, label=F, alpha=150, kind="se", conf=0.95)
legend("topright", legend=MDS.levels, pch=NMDS.symbols, col="black", inset=c(0.18, 0.0), cex=0.9, pt.bg=MDS.colors, pt.cex=1, bty="n", x.intersp = 1.2, y.intersp = 1)
par.new=T
plot(vars, p.max=0.05, cex=0.8, lwd=1)
##### save the figure and export to directory? ####
#dev.copy(pdf, "figures/Fig3.NMDS_MC.erratum.pdf", height=7, width=7, useDingbats=FALSE)
#dev.off()
##### NMDS plots for Porites compressa
##### ##### ##### ##### ##### ##### #####
##### ##### ##### ##### ##### #####
df.manova2.PC<-df.manova[(df.manova$Species=="PC" & df.manova$Period=="Bleaching"),]
NMDS<-subset(df.manova2.PC, select = c(-Time.point, -Period, -Site, -Species, -Condition, -Sample.ID, -Pair, -Depth.m))
colnames(NMDS)<-c("chl", "biomass", "proteins", "lipids", "carbs", "kJ", "d15N-H", "d13C-H", "C:N-H", "d15N-S", "d13C-S", "C:N-S")
NMDS.PC<-metaMDS(NMDS, distance='bray', k=2, trymax=100)
# stressplot(NMDS.PC)
# plot(NMDS.PC, type='t')
species.scores <- as.data.frame(scores(NMDS.PC, "species")) # extracts 'species' scores
species.scores$species <- rownames(species.scores) #creates a 'species columns' = this is the variables
factor.df<-df.manova2.PC[, c(2:8)]; names(factor.df)
NMDS.df<-data.frame(x=NMDS.PC$point[,1], y=NMDS.PC$point[,2],
Period=as.factor(factor.df[,1]),
Site=as.factor(factor.df[,2]),
Species=as.factor(factor.df[,3]),
Condition=as.factor(factor.df[,4]),
Sample.ID=as.factor(factor.df[,5]))
NMDS.df<-NMDS.df[, c(3:7,1,2)] # reorders columns
colnames(NMDS.df)[6:7]<-c("NMDS1", "NMDS2")
MDS.levels<-c("R44-B", "R44-NB", "R25-B", "R25-NB", "HIMB-B", "HIMB-NB")
MDS.colors<-c("lightpink", "firebrick2", "slategray1", "dodgerblue2","khaki", "darkorange")
NMDS.mean=aggregate(NMDS.df[,6:7],list(group=NMDS.df$Site:NMDS.df$Condition), mean)
NMDS.symbols<-c(25, 25, 23, 23, 21, 21)
groups<-NMDS.df$Site:NMDS.df$Condition
# vectors for variables
vars<-envfit(NMDS.PC, NMDS, permu=999)
#scores(vars, "vectors")
###################
# PLOTS
### Porites 2014-- Site:Condition
par(mfcol=c(2,2), mar=c(4,4,1,1), pty="sq")
NMDS.plot<-ordiplot(NMDS.PC, type='n', main=substitute(paste(italic("Porites compressa"))), cex.main=1, display="sites", xlim=c(-0.12, 0.12), ylim=c(-0.12, 0.12), cex.lab=0.8, cex.axis=0.8, xaxt='n', xlab="")
axis(side = 1, labels = FALSE, tck = -0.03)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
text(-0.066, 0.114, "Oct '14: Bleaching", cex=1)
points(NMDS.plot, "sites", pch=NMDS.symbols[groups], cex=1, col="black", bg=MDS.colors[groups])
with(NMDS.plot, points(NMDS.mean[,2], NMDS.mean[,3], pch=4))
ordiellipse(NMDS.plot, groups=NMDS.df$Site:NMDS.df$Condition, draw="polygon", col=MDS.colors, label=F, alpha=150, kind="se", conf=0.95)
##### ##### ##### ##### ##### ##### #####
##### ##### ##### ##### ##### #####
### Porites 2014-- Condition
MDS.levels<-c("Bleached", "Nonbleached")
MDS.colors<-c("lightpink", "slategray1")
NMDS.mean=aggregate(NMDS.df[,6:7],list(group=NMDS.df$Condition), mean)
NMDS.symbols<-c(21, 23)
groups2<-NMDS.df$Condition
NMDS.plot<-ordiplot(NMDS.PC, type='n', display="sites", xlim=c(-0.12, 0.12), ylim=c(-0.12, 0.12), cex.lab=0.8, cex.axis=0.8, ylab= "NMDS2", xlab="NMDS1")
axis(side = 1, labels = FALSE, tck = -0.03)
axis(side = 2, labels = FALSE, tck = -0.03)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
points(NMDS.plot, "sites", pch=NMDS.symbols[groups2], cex=1, col="black", bg=MDS.colors[groups2])
with(NMDS.plot, points(NMDS.mean[,2], NMDS.mean[,3], pch=4))
ordiellipse(NMDS.plot, groups=NMDS.df$Condition, draw="polygon", col=MDS.colors, label=F, alpha=150, kind="se", conf=0.95)
par.new=T
plot(vars, p.max=0.05, cex=0.8, lwd=1)
##### ##### ##### ##### ##### ##### #####
##### ##### ##### ##### ##### #####
### Porites 2015 ###
df.manova2.PC<-df.manova[(df.manova$Species=="PC" & df.manova$Period=="Recovery"),]
NMDS<-subset(df.manova2.PC, select = c(-Time.point, -Period, -Site, -Species, -Condition, -Sample.ID, -Pair, -Depth.m))
colnames(NMDS)<-c("chl", "biomass", "protein", "lipids", "carbs", "kJ", "d15N-host", "d13C-host", "C:N-host", "d15N-symb", "d13C-symb", "C:N-symb")
NMDS.PC<-metaMDS(NMDS, distance='bray', k=2, trymax=100)
# stressplot(NMDS.PC)
# plot(NMDS.PC, type='t')
species.scores <- as.data.frame(scores(NMDS.PC, "species")) # extracts 'species' scores
species.scores$species <- rownames(species.scores) #creates a 'species columns' = this is the variables
factor.df<-df.manova2.PC[, c(2:8)]; names(factor.df)
NMDS.df<-data.frame(x=NMDS.PC$point[,1], y=NMDS.PC$point[,2],
Period=as.factor(factor.df[,1]),
Site=as.factor(factor.df[,2]),
Species=as.factor(factor.df[,3]),
Condition=as.factor(factor.df[,4]),
Sample.ID=as.factor(factor.df[,5]))
NMDS.df<-NMDS.df[, c(3:7,1,2)] # reorders columns
colnames(NMDS.df)[6:7]<-c("NMDS1", "NMDS2")
MDS.levels<-c("R44-B", "R44-NB", "R25-B", "R25-NB", "HIMB-B", "HIMB-NB")
MDS.colors<-c("lightpink", "firebrick2", "slategray1", "dodgerblue2","khaki", "darkorange")
NMDS.mean=aggregate(NMDS.df[,6:7],list(group=NMDS.df$Site:NMDS.df$Condition), mean)
NMDS.symbols<-c(25, 25, 23, 23, 21, 21)
# vectors for variables
vars<-envfit(NMDS.PC, NMDS, permu=999)
# scores(vars, "vectors")
##### ##### ##### ##### ##### ##### #####
##### ##### ##### ##### ##### #####
### "Porites 2015---Site:Condition"
NMDS.plot<-ordiplot(NMDS.PC, type='n', display="sites", xlim=c(-0.12, 0.12), ylim=c(-0.12, 0.12), cex.lab=0.8, cex.axis=0.8, ylab="", xlab="", xaxt="n", yaxt="n")
axis(side = 1, labels = FALSE, tck = -0.03)
axis(side = 2, labels = FALSE, tck = -0.03)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
text(-0.066, 0.114, "Jan '15: Recovery", cex=1)
points(NMDS.plot, "sites", pch=NMDS.symbols[groups], cex=1, col="black", bg=MDS.colors[groups])
with(NMDS.plot, points(NMDS.mean[,2], NMDS.mean[,3], pch=4))
ordiellipse(NMDS.plot, groups=NMDS.df$Site:NMDS.df$Condition, draw="polygon", col=MDS.colors, label=F, alpha=150, kind="se", conf=0.95)
legend("topright", legend=MDS.levels, pch=NMDS.symbols, col="black", cex=0.85, inset=c(0.06,0), pt.bg=MDS.colors, pt.cex=1, bty="n", x.intersp = 1.2, y.intersp = 1)
##### ##### ##### ##### ##### ##### #####
##### ##### ##### ##### ##### #####
### "Porites 2015---Condition"
MDS.levels<-c("Bleached", "Non-bleached")
MDS.colors<-c("lightpink", "slategray1")
NMDS.mean=aggregate(NMDS.df[,6:7],list(group=NMDS.df$Condition), mean)
NMDS.symbols<-c(21, 23)
NMDS.plot<-ordiplot(NMDS.PC, type='n', display="sites", xlim=c(-0.12, 0.12), ylim=c(-0.12, 0.12), cex.lab=0.8, cex.axis=0.8, ylab="", yaxt='n')
axis(side = 2, labels = FALSE, tck = -0.03)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
#text(-0.15, 0.11, "d", cex=1.5)
points(NMDS.plot, "sites", pch=NMDS.symbols[groups2], cex=1, col="black", bg=MDS.colors[groups2])
with(NMDS.plot, points(NMDS.mean[,2], NMDS.mean[,3], pch=4))
ordiellipse(NMDS.plot, groups=NMDS.df$Condition, draw="polygon", col=MDS.colors, label=F, alpha=150, kind="se", conf=0.95)
legend("topright", legend=MDS.levels, pch=NMDS.symbols, col="black", inset=c(0.18, 0.0), cex=0.9, pt.bg=MDS.colors, pt.cex=1, bty="n", x.intersp = 1.2, y.intersp = 1)
par.new=T
plot(vars, p.max=0.05, cex=0.8, lwd=1)
##### save the figure and export to directory? ####
#dev.copy(pdf, "figures/Fig4.NMDS_PC.erratum.pdf", height=7, width=7, useDingbats=FALSE)
#dev.off()
##### MANOVA
MANOVA outputs These outputs from MANOVA cane be found in Table S3.
The first MANOVA shows Montipora capitata
# Species sufficently distinct (NMDS shows this), run MANOVA with Species separated
# from above---df.manova2.MC is the df with Montipora alone
MC.Manova<-df.manova[(df.manova$Species=="MC"),]
MC.Manova.data<-MC.Manova[, c(9:20)] # makes df without factor columns
# from above---df.manova2.PC is the df with Porites alone
PC.Manova<-df.manova[(df.manova$Species=="PC"),]
PC.Manova.data<-PC.Manova[, c(9:20)] # makes df without factor columns
#### #### #### #### #### ####
#### #### #### #### #### ####
#### Montipora
MVA.MC<-adonis2(MC.Manova.data~Period*Site*Condition, data=MC.Manova,permutations=1000,
method="bray", sqrt.dist = TRUE)
MVA.MC # effect of Period > Condition > Site and Period:Condition## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = MC.Manova.data ~ Period * Site * Condition, data = MC.Manova, permutations = 1000, method = "bray", sqrt.dist = TRUE)
## Df SumOfSqs R2 F Pr(>F)
## Period 1 0.5611 0.16053 11.7135 0.000999 ***
## Site 2 0.2207 0.06313 2.3031 0.011988 *
## Condition 1 0.1491 0.04266 3.1126 0.008991 **
## Period:Site 2 0.0881 0.02521 0.9198 0.497502
## Period:Condition 1 0.1088 0.03112 2.2708 0.043956 *
## Site:Condition 2 0.1016 0.02906 1.0601 0.343656
## Period:Site:Condition 2 0.0625 0.01788 0.6522 0.878122
## Residual 46 2.2037 0.63042
## Total 57 3.4956 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# posthocs
MC.bray=vegdist(MC.Manova.data) #matrix as a distance matrix
# pairwise.perm.manova(MC.bray, MC.Manova$Period, nperm=1000, p.method="holm") # Periods differ (p=0.001)
# pairwise.perm.manova(MC.bray, MC.Manova$Site, nperm=1000, p.method="holm") # Sites not different (p>0.11)
# pairwise.perm.manova(MC.bray, MC.Manova$Condition, nperm=1000, p.method="holm") # not different (p=0.083)
# pairwise.perm.manova(MC.bray, MC.Manova$Period:MC.Manova$Condition, nperm=1000, p.method="holm")
# @ bleaching--bleached corals different from non-bleached
# @ recovery--bleached corals not different fron non-bleached (p=0.305)
#### #### #### #### #### ####
#### #### #### #### #### #### The second MANOVA shows Porites compressa
#### Porites
MVA.PC<-adonis2(PC.Manova.data~Period*Site*Condition, data=PC.Manova, permutations=1000,
method="bray", sqrt.dist = TRUE)
MVA.PC # effect of Period > Condition > Period:Condition## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = PC.Manova.data ~ Period * Site * Condition, data = PC.Manova, permutations = 1000, method = "bray", sqrt.dist = TRUE)
## Df SumOfSqs R2 F Pr(>F)
## Period 1 0.4971 0.13541 10.0244 0.000999 ***
## Site 2 0.1414 0.03852 1.4257 0.110889
## Condition 1 0.1904 0.05187 3.8399 0.000999 ***
## Period:Site 2 0.1448 0.03943 1.4594 0.100899
## Period:Condition 1 0.1430 0.03894 2.8825 0.007992 **
## Site:Condition 2 0.1333 0.03631 1.3438 0.134865
## Period:Site:Condition 2 0.0904 0.02463 0.9116 0.530470
## Residual 47 2.3309 0.63490
## Total 58 3.6713 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# posthocs
PC.bray=bcdist(PC.Manova.data) #matrix as a distance matrix
# pairwise.perm.manova(PC.bray, PC.Manova$Period, nperm=1000, p.method="holm") # different (p=0.001)
# pairwise.perm.manova(PC.bray, PC.Manova$Condition, nperm=1000, p.method="holm") # different (p=0.005)
# pairwise.perm.manova(PC.bray, PC.Manova$Period:PC.Manova$Condition, nperm=1000, p.method="holm")
# @ bleaching--bleached corals different from non-bleached
# @ recovery--bleached corals not different fron non-bleached (p=0.110)Univariate analysis
separate models by species Here, each response variable is now analyzed in a linear mixed effect model. For each response variable model selection is performed as follows:
Divide the species The NMDS showed species were distinct, dividing the species can allow for univariate tests to look at the species-specifc responses more clearly.
Using LMER First, check for valid assumptions of ANVOA and normal distribution using graphical inspections of residuals and QQ plots.
Linear mixed effect models are run with all fixed effects (Period, Site, Condition) and the repeated sampling of individuals is accounted for using coral ID (Sample ID) and colony pairs (Pairs) as a random effects. This model is labeled as: full.test. This model is fit with restricted maximum likelihood (REML). Inspect the output and identify interactive effects that can be dropped. Next, compare models with different fixed effect.
Model Selection After fitting the first model (full), compare to a reduced model where non-significant interactive terms have been dropped. This model is labeled: add. To compare models, the full model must be refit using maximum likelihood, set by REML=FALSE, in order to compare models with different fixed effects. This re-fit ML-model– for comparison only– is called: full. With both models fit by maximum likelihood, use likelihood ratio tests and compare the two models and see if there is reduced model (e.g., add) differs from full model. If not (P > 0.15, since chi-square is usually inflated), drop non-significant effects and proceed with the best model.
Fit Final Model Run the best model with end tag xxx.mod.FINAL, fit with REML, generate ANVOA table (Type II SS), and run posthoc test (emmeans) where necessary. Posthoc tests should be run as ‘slices’, where a comparison is made among groups and primary effects tested.
Output from these models can be found summarized in Tables 1, S4 (MC physiology), S5 (PC physiology).
Montipora capitata
Code will test for assumptions of ANOVA. Where assumptions are not met, data transformations are applied, with help selecting an appropriate transformation aided with a Box-Cox test. A loop generates diagnostic plots (QQ plots, histogram of residuals) and tests metrics for adherence to ANOVA assumptions.
Linear Mixed Effect Models
Total chlorophyll (ug cm-2)
###################
# chl a+c2
###################
full.test<-lmer(chl~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.MC, na.action=na.exclude)
anova(full.test, type=2) # only slight effect of Site
#### compare models ####
# full model
full<-lmer(chl~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.MC, na.action=na.exclude)
# additive model
add<-lmer(chl~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.MC, na.action=na.exclude)
anova(full, add) # no difference# final model
Chl.mod.FINAL<-lmer(chl~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.MC, na.action=na.exclude)
# ANOVA table
anova(Chl.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 235.380 235.380 1 28 65.2446 8.541e-09 ***
## Site 12.354 6.177 2 12 1.7122 0.221733
## Condition 54.080 54.080 1 14 14.9905 0.001693 **
## Period:Condition 16.541 16.541 1 28 4.5851 0.041097 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(Chl.mod.FINAL, pairwise~Condition|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Condition emmean SE df lower.CL upper.CL .group
## B 1.95 0.594 45.8 0.757 3.15 a
## NB 5.41 0.594 45.8 4.219 6.61 b
##
## Period = Recovery:
## Condition emmean SE df lower.CL upper.CL .group
## B 6.96 0.594 45.8 5.768 8.16 a
## NB 8.32 0.594 45.8 7.130 9.52 a
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(Chl.mod.FINAL), ylab="total chlorophyll ug/cm2")
# Bleaching: B and NB differ (a vs. b)
# Recovery: B and NB do not differ (a vs. a)Total biomass (AFDW mg cm-2)
###################
# AFDW biomass
###################
full.test<-lmer(biomass~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.MC, na.action=na.exclude)
anova(full.test, type=2) # only slight effect of Site
#### compare models ####
# full model
full<-lmer(biomass~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.MC, na.action=na.exclude)
# additive model
add<-lmer(biomass~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.MC, na.action=na.exclude)
anova(full, add) # no difference# final model
AFDW.mod.FINAL<-lmer(biomass~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.MC, na.action=na.exclude)
# rand(AFDW.mod.FINAL)
# ANOVA table
print(anova(AFDW.mod.FINAL, type=2), digits=6)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 2912.240 2912.240 1 28 55.48377 4.111e-08 ***
## Site 145.006 72.503 2 26 1.38132 0.269079
## Condition 103.915 103.915 1 26 1.97978 0.171256
## Period:Condition 394.610 394.610 1 28 7.51808 0.010524 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(AFDW.mod.FINAL, pairwise~Condition|Period)
CLD(posthoc, Letters=letters); ## Period = Bleaching:
## Condition emmean SE df lower.CL upper.CL .group
## B 19.2 2.06 51.5 15.0 23.3 a
## NB 27.4 2.06 51.5 23.3 31.5 b
##
## Period = Recovery:
## Condition emmean SE df lower.CL upper.CL .group
## NB 36.2 2.06 51.5 32.1 40.4 a
## B 38.2 2.06 51.5 34.1 42.3 a
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05# Bleaching period: B < NB (a vs. b)
# Recovery period: B = NB (a vs a)
plot(allEffects(AFDW.mod.FINAL), ylab="M. capitata AFDW mg/cm2")
Protein (g gdw-1)
###################
# prot.gdw
###################
full.test<-lmer(prot~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.MC, na.action=na.exclude)
anova(full.test, type=2) # only effect of Period
#### compare models ####
# full model
full<-lmer(prot~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.MC, na.action=na.exclude)
# additive model
add<-lmer(prot~Period+Site+Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.MC, na.action=na.exclude)
anova(full, add) # no difference# final model
Prot.gdw.mod.FINAL<-lmer(prot~Period+Site+Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.MC, na.action=na.exclude)
# ANOVA table
anova(Prot.gdw.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 0.0031835 0.0031835 1 28.362 7.6886 0.009713 **
## Site 0.0006606 0.0003303 2 25.972 0.7978 0.461060
## Condition 0.0005335 0.0005335 1 26.001 1.2886 0.266662
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(Prot.gdw.mod.FINAL)
posthoc<-emmeans(Prot.gdw.mod.FINAL, pairwise~Period)
CLD(posthoc, Letters=letters);## Period emmean SE df lower.CL upper.CL .group
## Recovery 0.0601 0.00443 26 0.0510 0.0693 a
## Bleaching 0.0750 0.00443 26 0.0659 0.0842 b
##
## Results are averaged over the levels of: Site, Condition
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(Prot.gdw.mod.FINAL), ylab="M. capitata protein/gdw")
Lipids (mg gdw-1)
###################
# lip.gdw
###################
full.test<-lmer(lipid~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.MC, na.action=na.exclude)
anova(full.test, type=2) # only slight effect of Site
#### compare models ####
# full model
full<-lmer(lipid~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.MC, na.action=na.exclude)
# additive model
add<-lmer(lipid~Period+Site+Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.MC, na.action=na.exclude)
anova(full, add) # no difference# final model
Lipid.gdw.mod.FINAL<-lmer(lipid~Period+Site+Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.MC, na.action=na.exclude)
# ANOVA table
anova(Lipid.gdw.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 0.0006620 0.0006620 1 27.114 0.1628 0.68974
## Site 0.0212351 0.0106175 2 25.014 2.6113 0.09333 .
## Condition 0.0047278 0.0047278 1 25.027 1.1628 0.29117
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(Lipid.gdw.mod.FINAL)
plot(allEffects(Lipid.gdw.mod.FINAL), ylab="M. capitata lipids/gdw")
Carbohydrates (g gdw-1)
###################
# carb.gdw
###################
# full model
full.test<-lmer(carb~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.MC, na.action=na.exclude)
anova(full.test, type=2) # only slight effect of Period
#### compare models ####
# full model
full<-lmer(carb~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.MC, na.action=na.exclude)
# additive model
add<-lmer(carb~Period+Site+Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.MC, na.action=na.exclude)
anova(full, add) # no difference# final model
Carbs.gdw.mod.FINAL<-lmer(carb~Period+Site+Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.MC, na.action=na.exclude)
# ANOVA table
anova(Carbs.gdw.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 0.0028233 0.00282327 1 28.107 3.4811 0.07254 .
## Site 0.0032777 0.00163885 2 26.097 2.0207 0.15279
## Condition 0.0001583 0.00015827 1 26.106 0.1951 0.66230
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(Carbs.gdw.mod.FINAL)
plot(allEffects(Carbs.gdw.mod.FINAL), ylab="M. capitata carbs/gdw")
Biomass Kilojoule (kJ gdw-1)
###################
# kJ.gdw
###################
full.test<-lmer(kJ~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.MC, na.action=na.exclude)
anova(full.test, type=2) # only slight effect of Site
#### compare models ####
# full model
full<-lmer(kJ~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.MC, na.action=na.exclude)
# additive model
add<-lmer(kJ~Period+Site+Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.MC, na.action=na.exclude)
anova(full, add) # no difference# final model
kJ.gdw.mod.FINAL<-lmer(kJ~Period+Site+Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.MC, na.action=na.exclude)
# ANOVA table
anova(kJ.gdw.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 0.039 0.0389 1 27.369 0.0063 0.93740
## Site 37.867 18.9336 2 25.082 3.0595 0.06472 .
## Condition 5.666 5.6662 1 25.106 0.9156 0.34776
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(kJ.gdw.mod.FINAL)
plot(allEffects(kJ.gdw.mod.FINAL), ylab="M. capitata kJ/gdw")
Porites compressa
Linear Mixed Effect Models
Chlorophylls (ug cm-2)
###################
#chl a+c2
################### SITE RETAINED-- Period, Condition, Period:Condition effect, almost a Site effect
# full model
full<-lmer(chl~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.PC, na.action=na.exclude)
anova(full, type=2) # lots of effects, use full model# final model
PC.tot.chl.mod.FINAL<-lmer(chl~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.PC, na.action=na.exclude)
# ANOVA table
print(anova(PC.tot.chl.mod.FINAL, type=2), digits =6)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 1641.570 1641.570 1 23.9998 258.38515 2.3854e-14 ***
## Site 61.579 30.789 2 12.0000 4.84628 0.028657 *
## Condition 257.603 257.603 1 12.0001 40.54699 3.5693e-05 ***
## Period:Site 39.845 19.923 2 23.9998 3.13586 0.061673 .
## Period:Condition 187.462 187.462 1 23.9998 29.50682 1.3979e-05 ***
## Site:Condition 86.013 43.007 2 12.0001 6.76929 0.010762 *
## Period:Site:Condition 8.942 4.471 2 23.9998 0.70376 0.504647
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(PC.tot.chl.mod.FINAL)
posthoc<-emmeans(PC.tot.chl.mod.FINAL, pairwise~Condition|Period:Site)
CLD(posthoc, Letters=letters)## Period = Bleaching, Site = HIMB:
## Condition emmean SE df lower.CL upper.CL .group
## B 1.01 1.28 45.2 -1.56 3.59 a
## NB 12.99 1.28 45.2 10.42 15.57 b
##
## Period = Bleaching, Site = Reef 25:
## Condition emmean SE df lower.CL upper.CL .group
## B 1.76 1.28 45.2 -0.81 4.34 a
## NB 8.16 1.28 45.2 5.58 10.73 b
##
## Period = Bleaching, Site = Reef 44:
## Condition emmean SE df lower.CL upper.CL .group
## B 1.29 1.28 45.2 -1.28 3.87 a
## NB 8.08 1.28 45.2 5.51 10.65 b
##
## Period = Recovery, Site = HIMB:
## Condition emmean SE df lower.CL upper.CL .group
## B 15.23 1.28 45.2 12.66 17.80 a
## NB 20.43 1.28 45.2 17.86 23.01 b
##
## Period = Recovery, Site = Reef 25:
## Condition emmean SE df lower.CL upper.CL .group
## NB 11.92 1.28 45.2 9.34 14.49 a
## B 14.62 1.28 45.2 12.04 17.19 a
##
## Period = Recovery, Site = Reef 44:
## Condition emmean SE df lower.CL upper.CL .group
## B 16.21 1.28 45.2 13.64 18.79 a
## NB 17.66 1.28 45.2 15.08 20.23 a
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05posthoc<-emmeans(PC.tot.chl.mod.FINAL, pairwise~Condition|Site)
CLD(posthoc, Letters=letters)## Site = HIMB:
## Condition emmean SE df lower.CL upper.CL .group
## B 8.12 0.999 23.6 6.06 10.2 a
## NB 16.71 0.999 23.6 14.65 18.8 b
##
## Site = Reef 25:
## Condition emmean SE df lower.CL upper.CL .group
## B 8.19 0.999 23.6 6.13 10.3 a
## NB 10.04 0.999 23.6 7.97 12.1 a
##
## Site = Reef 44:
## Condition emmean SE df lower.CL upper.CL .group
## B 8.75 0.999 23.6 6.69 10.8 a
## NB 12.87 0.999 23.6 10.80 14.9 b
##
## Results are averaged over the levels of: Period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(PC.tot.chl.mod.FINAL), ylab="P. compressa chl ug/cm2")
Total biomass (AFDW mg cm-2)
###################
# AFDW biomass
###################
# full model
full.test<-lmer(biomass~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.PC, na.action=na.exclude)
anova(full, type=2) # only Condition effect, drop interactions
#### compare models ####
full<-lmer(biomass~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.PC, na.action=na.exclude)
# additive model
add<-lmer(biomass~Period+Site+Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.PC, na.action=na.exclude)
# compare models
AIC(full, add)
anova(full, add) # no difference in effects# final model
PC.AFDW.mod<-lmer(biomass~Period+Site+Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.PC, na.action=na.exclude)
# ANOVA table
print(anova(PC.AFDW.mod, type=2), digits=6)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 248.036 248.036 1 55 1.90991 0.172561
## Site 74.351 37.176 2 55 0.28626 0.752182
## Condition 691.067 691.067 1 55 5.32130 0.024858 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(PC.AFDW.mod)
posthoc<-emmeans(PC.AFDW.mod, pairwise~Condition|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Condition emmean SE df lower.CL upper.CL .group
## B 33.8 2.55 47.6 28.7 39.0 a
## NB 40.6 2.55 47.6 35.5 45.8 b
##
## Period = Recovery:
## Condition emmean SE df lower.CL upper.CL .group
## B 37.9 2.55 47.6 32.8 43.0 a
## NB 44.7 2.55 47.6 39.6 49.8 b
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05# @ Bleaching: B < NB, but @ Recovery B = NB
plot(allEffects(PC.AFDW.mod), ylab="P. compressa AFDW mg/cm2")
Protein (g gdw-1)
###################
# prot.gdw
###################
full.test<-lmer(prot~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.PC, na.action=na.exclude)
anova(full.test, type=2) # only Period:Condition effect
#### compare models ####
# full model
full<-lmer(prot~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.PC, na.action=na.exclude)
# additive model
add<-lmer(prot~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.PC, na.action=na.exclude)
anova(full, add) # no difference# final model
PC.prot.mod.FINAL<-lmer(prot~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.PC, na.action=na.exclude)
# ANOVA table
anova(PC.prot.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 0.0000891 0.0000891 1 27.649 0.0856 0.77203
## Site 0.0029600 0.0014800 2 25.890 1.4223 0.25939
## Condition 0.0000761 0.0000761 1 25.892 0.0731 0.78904
## Period:Condition 0.0076778 0.0076778 1 27.666 7.3784 0.01125 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(PC.prot.mod.FINAL)
posthoc<-emmeans(PC.prot.mod.FINAL, pairwise~Condition|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Condition emmean SE df lower.CL upper.CL .group
## NB 0.107 0.00925 50.2 0.0884 0.126 a
## B 0.128 0.00962 50.8 0.1082 0.147 a
##
## Period = Recovery:
## Condition emmean SE df lower.CL upper.CL .group
## B 0.102 0.00925 50.2 0.0832 0.120 a
## NB 0.127 0.00925 50.2 0.1084 0.146 a
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(PC.prot.mod.FINAL), ylab="P. compressa protein/gdw")
Lipids (g gdw-1)
###################
# lip.gdw
###################
full.test<-lmer(lipid~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.PC, na.action=na.exclude)
anova(full, type=2) # only Period, Period:Site effects, drop interactions
#### compare models ####
# full model
full<-lmer(lipid~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.PC, na.action=na.exclude)
# additive model
add<-lmer(lipid~Period+Site+Condition+Period:Site+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.PC, na.action=na.exclude)
anova(full, add) # no difference# final model
PC.lipids.mod.FINAL<-lmer(lipid~Period+Site+Condition+Period:Site+(1|Sample.ID)+(1|Pair), data=phys.datatests.PC, na.action=na.exclude)
# ANOVA table
anova(PC.lipids.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 0.064354 0.064354 1 52 12.184 0.0009915 ***
## Site 0.004764 0.002382 2 52 0.451 0.6394621
## Condition 0.009803 0.009803 1 52 1.856 0.1789643
## Period:Site 0.059370 0.029685 2 52 5.620 0.0061702 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(PC.lipids.mod.FINAL)
posthoc<-emmeans(PC.lipids.mod.FINAL, pairwise~Site|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Site emmean SE df lower.CL upper.CL .group
## HIMB 0.386 0.0230 32.4 0.339 0.432 a
## Reef 44 0.431 0.0230 32.4 0.384 0.478 a
## Reef 25 0.440 0.0244 34.9 0.390 0.489 a
##
## Period = Recovery:
## Site emmean SE df lower.CL upper.CL .group
## Reef 44 0.320 0.0230 32.4 0.273 0.367 a
## Reef 25 0.327 0.0230 32.4 0.281 0.374 a
## HIMB 0.408 0.0230 32.4 0.361 0.455 b
##
## Results are averaged over the levels of: Condition
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05plot(allEffects(PC.lipids.mod.FINAL), ylab="P. compressa lipids/gdw")
Carbohydrates (g gdw-1)
###################
# carb.gdw
###################
full.test<-lmer(carb~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.PC, na.action=na.exclude)
anova(full.test) # no effects
#### compare models ####
# full model
full<-lmer(carb~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.PC, na.action=na.exclude)
# additive model
add<-lmer(carb~Period+Site+Condition+Period:Site+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.PC, na.action=na.exclude)
anova(full, add) # no difference# final model
PC.carb.mod.FINAL<-lmer(carb~Period+Site+Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.PC, na.action=na.exclude)
# ANOVA table
anova(PC.carb.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 0.00141663 0.00141663 1 28.838 2.6527 0.1143
## Site 0.00063323 0.00031661 2 26.075 0.5929 0.5600
## Condition 0.00125195 0.00125195 1 26.082 2.3443 0.1378# rand(PC.carb.mod.FINAL)
plot(allEffects(PC.carb.mod.FINAL), ylab="P. compressa carbs/gdw")
Biomass kilojoule (kJ gdw-1)
###################
# kJ.gdw
###################
full.test<-lmer(kJ~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=phys.datatests.PC, na.action=na.exclude)
anova(full, type=2) # only Period, Period:Site effects, drop interactions
#### compare models ####
# full model
full<-lmer(kJ~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.PC, na.action=na.exclude)
# additive model
add<-lmer(kJ~Period+Site+Condition+Period:Site+(1|Sample.ID)+(1|Pair), REML=FALSE, data=phys.datatests.PC, na.action=na.exclude)
anova(full, add) # no difference# final model
PC.kJ.mod.FINAL<-lmer(kJ~Period+Site+Condition+Period:Site+(1|Sample.ID)+(1|Pair), data=phys.datatests.PC, na.action=na.exclude)
# ANOVA table
anova(PC.kJ.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 92.140 92.140 1 52 14.0711 0.0004442 ***
## Site 4.004 2.002 2 52 0.3057 0.7378908
## Condition 8.596 8.596 1 52 1.3127 0.2571540
## Period:Site 69.078 34.539 2 52 5.2746 0.0082094 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(PC.kJ.mod.FINAL)
posthoc<-emmeans(PC.kJ.mod.FINAL, pairwise~Site|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Site emmean SE df lower.CL upper.CL .group
## HIMB 19.4 0.809 32.4 17.7 21.0 a
## Reef 44 20.9 0.809 32.4 19.2 22.5 a
## Reef 25 21.8 0.861 34.9 20.0 23.5 a
##
## Period = Recovery:
## Site emmean SE df lower.CL upper.CL .group
## Reef 44 17.2 0.809 32.4 15.5 18.8 a
## Reef 25 17.4 0.809 32.4 15.7 19.0 a
## HIMB 19.9 0.809 32.4 18.2 21.5 a
##
## Results are averaged over the levels of: Condition
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05plot(allEffects(PC.kJ.mod.FINAL), ylab="P. compressa kJ/gdw")
Physiology Figures
Figures shown here are total chlorophyll, total biomass and found in Figure 5
# below is for x-axis labels
Period.Reef<-paste(expression("R44\nBleaching", "R25\nBleaching", "HIMB\nBleaching", "R44\nRecovery", "R25\nRecovery", "HIMB\nRecovery"))
df.MC$Site<-factor(df.MC$Site, levels=c("Reef 44", "Reef 25", "HIMB")) # for ordering x axis
# figure formatting script
Fig.formatting<-(theme_classic()) +
theme(text=element_text(size=8),
axis.line=element_blank(),
legend.justification=c(1,1), legend.position=c(1,1),
legend.background = element_rect("transparent", colour=NA),
legend.text.align = 0,
legend.text=element_text(size=8),
legend.title = element_blank(),
panel.border = element_rect(fill=NA, colour = "black", size=1),
aspect.ratio=0.7,
axis.ticks.length=unit(0.25, "cm"),
axis.text.y=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10),
axis.text.x=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=6)) +
theme(legend.key.size = unit(0.4, "cm")) +
theme(panel.spacing=unit(c(0, 0, 0, 0), "cm"))
pd <- position_dodge(0.71) #offset for error bars and columns
###############################
# Montipora capitata figures #
###############################
########## chlorophyll (a + c2) #############
Fig.MC.chl<-ggplot(df.MC, aes(x=Period:Site, y=chl, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=chl-chlSE, ymax=chl+chlSE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("Chlorophyll", ~italic("a"+"c"[2]), ~(mu*g~cm^-2), sep=""))) +
scale_y_continuous(expand=c(0,0), limits=c(0, 20)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) +
geom_vline(xintercept = 3.5, linetype="dotted")+ Fig.formatting
########## Tissue biomass ##########
Fig.MC.AFDW <-ggplot(df.MC, aes(x=Period:Site, y=AFDW, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=AFDW-AFDWSE, ymax=AFDW+AFDWSE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("Tissue biomass", ~(mg~cm^-2), sep=""))) +
scale_y_continuous(expand=c(0,0), limits=c(0, 80)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached)))+ guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted")+
Fig.formatting
########## Proteins gdw ##########
Fig.MC.prot.gdw<-ggplot(df.MC, aes(x=Period:Site, y=prot, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=prot-protSE, ymax=prot+protSE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("Proteins", ~(g~gdw^-1), sep=""))) +
scale_y_continuous(expand=c(0,0), limits=c(0, 0.2)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) +
geom_vline(xintercept = 3.5, linetype="dotted")+ Fig.formatting
########## Carbohydrates gdw ##########
Fig.MC.carbs<-ggplot(df.MC, aes(x=Period:Site, y=carbs, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=carbs-carbsSE, ymax=carbs+carbsSE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("Carbohydrates", ~(g~gdw^-1), sep=""))) +
scale_y_continuous(expand=c(0,0), limits=c(0, 0.2)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached)))+ guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted")+ Fig.formatting
########## Lipids gdw ##########
Fig.MC.lipid<-ggplot(df.MC, aes(x=Period:Site, y=lipids, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=lipids-lipidsSE, ymax=lipids+lipidsSE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("Lipids", ~(g~gdw^-1), sep=""))) +
scale_y_continuous(expand=c(0,0), limits=c(0, 0.8)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached)))+ guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted")+ theme(legend.position="none") +
Fig.formatting
########## Energy gdw ##########
Fig.MC.kJ<-ggplot(df.MC, aes(x=Period:Site, y=kJ, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=kJ-kJSE, ymax=kJ+kJSE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("Energy", ~(kJ~gdw^-1), sep=""))) +
scale_y_continuous(expand=c(0,0), limits=c(0, 40)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached)))+ guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted")+ theme(legend.position="none") +
Fig.formatting
##########
##########
df.PC$Site<-factor(df.PC$Site, levels=c("Reef 44", "Reef 25", "HIMB"))
###############################
# Porites compressa figures #
###############################
########## chlorophyll (a + c2) #############
Fig.PC.chl<-ggplot(df.PC, aes(x=Period:Site, y=chl, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=chl-chlSE, ymax=chl+chlSE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("Chlorophyll", ~italic("a"+"c"[2]), ~(mu*g~cm^-2), sep=""))) +
scale_y_continuous(expand=c(0,0), limits=c(0, 30)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) + guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted") + Fig.formatting
########## Tissue biomass ##########
Fig.PC.AFDW <-ggplot(df.PC, aes(x=Period:Site, y=AFDW, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=AFDW-AFDWSE, ymax=AFDW+AFDWSE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("Tissue biomass", ~(mg~cm^-2), sep=""))) +
scale_y_continuous(expand=c(0,0), limits=c(0, 80)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) + guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted")+
Fig.formatting
########## Proteins gdw ##########
Fig.PC.prot.gdw<-ggplot(df.PC, aes(x=Period:Site, y=prot, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=prot-protSE, ymax=prot+protSE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("Proteins", ~(g~gdw^-1), sep=""))) +
scale_y_continuous(expand=c(0,0), limits=c(0, 0.2)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) + guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted") +
Fig.formatting
########## Carbohydrates gdw ##########
Fig.PC.carbs<-ggplot(df.PC, aes(x=Period:Site, y=carbs, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=carbs-carbsSE, ymax=carbs+carbsSE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("Carbohydrates", ~(g~gdw^-1), sep=""))) +
scale_y_continuous(expand=c(0,0), limits=c(0, 0.2)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached)))+ guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted")+
Fig.formatting
########## Lipids gdw ##########
Fig.PC.lipid<-ggplot(df.PC, aes(x=Period:Site, y=lipids, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=lipids-lipidsSE, ymax=lipids+lipidsSE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("Lipids", ~(g~gdw^-1), sep=""))) +
scale_y_continuous(expand=c(0,0), limits=c(0, 0.8)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached)))+ guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted")+
Fig.formatting
########## Energy gdw ##########
Fig.PC.kJ<-ggplot(df.PC, aes(x=Period:Site, y=kJ, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=kJ-kJSE, ymax=kJ+kJSE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("Energy", ~(kJ~gdw^-1), sep=""))) +
scale_y_continuous(expand=c(0,0), limits=c(0, 40)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached)))+ guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted") +
Fig.formatting
plot_grid(Fig.MC.chl, Fig.PC.chl,
Fig.MC.AFDW, Fig.PC.AFDW,
labels = c("a", "b", "c", "d"), ncol = 2)
Figure 5. Physiology figures are mean +/-SE with Montipora capitata* on the left, Porites compressa on the right*.
Figures here are for biomass composition/energy reserves and are total proteins, lipids, carbohydrates, and energy content. These figures can be found in manuscript Figure 6.
plot_grid(Fig.MC.prot.gdw, Fig.PC.prot.gdw,
Fig.MC.lipid, Fig.PC.lipid,
labels = c("a", "b", "c", "d"), ncol = 2
)
Figure 6. Physiology figures are mean +/-SE with Montipora capitata* on the left, Porites compressa on the right*
plot_grid(Fig.MC.carbs, Fig.PC.carbs,
Fig.MC.kJ, Fig.PC.kJ,
labels = c("e", "f", "g", "h"), ncol = 2)
Figure 6. Physiology figures are mean +/-SE with Montipora capitata* on the left, Porites compressa on the right*
##Stable isotopes Isotope measurements are for isolated coral host and algal symbiont (Symbiodiniaceae) tissues. Carbon (δ13C) and nitrogen (δ15N) isotopic values and molar ratios of carbon:nitrogen (C:N) for coral host (δ13CH, δ15NH, C:NH) and algal symbiont (δ13CS, δ15NS, C:NS) tissues were determined with a Costech elemental combustion system coupled to a Thermo-Finnigan Delta Plus XP Isotope Ratio Mass-Spectrometer. Isotopic data are reported in delta values (δ) using the conventional permil (‰) notation and expressed relative to Vienna Pee-Dee Belemnite (V-PBD) and atmospheric N2 standards (air) for carbon and nitrogen, respectively.
Linear mixed effect models and anlaysis are separated by Species. First showing carbon, then nitrogen for animal and symbiont tissues in Montipora capitata followed by Porites compressa
Linear Mixed Effect Models
Montipora capitata
Results from these models are summarized in Table 1 and in Table S6
# iso.data -- all isotope data
# remove and order columns
iso.df<-iso.data[ , !names(iso.data) %in% c("host.ug.N", "host.ug.C", "symb.ug.N", "symb.ug.C")]
# MC data only
MC.iso.data<-iso.df[(iso.df$Species=="MC"),]
# data for transformations and analyses
MC.iso.tests<-MC.iso.data
## tests for MC isotopes data
for(i in c(10:17)){
Y<-MC.iso.tests[,i]
full<-lmer(Y~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
R <- resid(full) #save glm residuals
Y.shapiro <- shapiro.test(R) #runs a normality test on residuals
print(Y.shapiro) # null = normally distrubuted (P<0.05 = non-normal
op<-par(mfrow = c(2,2), mar=c(5,4,1,2), pty="sq")
plot(full, add.smooth = FALSE, which=1)
QQ <- qqnorm(R, main = colnames(MC.iso.tests)[i])
QQline <- qqline(R)
hist(R, xlab="Residuals", main = colnames(MC.iso.tests)[i])
}# boxcox(symb.d15N~Period*Site*Condition, data=MC.iso.tests, lambda=seq(-1, 2, length=10))
MC.iso.tests$host.d15N<-(MC.iso.tests$host.d15N)^0.25
MC.iso.tests$symb.C.N<-(MC.iso.tests$symb.C.N)^0.5
MC.iso.tests$symb.d15N<-(MC.iso.tests$symb.d15N)^0.25M. capitata, Carbon isotope, host: δ13CH
#############
# host d13C
#############
full.test<-lmer(host.d13C~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
anova(full, type=2) # only Site effect
#### compare models ####
# full model
full<-lmer(host.d13C~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(host.d13C~Period+Site+Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
MC.host.d13C.mod.FINAL<-lmer(host.d13C~Period+Site+Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(MC.host.d13C.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 0.6000 0.6000 1 29 1.1083 0.30115
## Site 3.5474 1.7737 2 12 3.2762 0.07323 .
## Condition 3.5613 3.5613 1 14 6.5782 0.02246 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(MC.host.d13C.mod.FINAL)
posthoc<-emmeans(MC.host.d13C.mod.FINAL, pairwise~Condition|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Condition emmean SE df lower.CL upper.CL .group
## NB -15.9 0.282 24.6 -16.5 -15.3 a
## B -15.2 0.282 24.6 -15.8 -14.6 b
##
## Period = Recovery:
## Condition emmean SE df lower.CL upper.CL .group
## NB -15.7 0.282 24.6 -16.3 -15.1 a
## B -15.0 0.282 24.6 -15.6 -14.4 b
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(MC.host.d13C.mod.FINAL), ylab="Mcap host d13C")
M. capitata, Carbon isotope, symbiont: δ13CS
#############
# symb d13C
#############
full.test<-lmer(symb.d13C~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
anova(full, type=2) # Period and Site effect
#### compare models ####
# full model
full<-lmer(symb.d13C~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(symb.d13C~Period+Site+Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
MC.symb.d13C.mod.FINAL<-lmer(symb.d13C~Period+Site+Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(MC.symb.d13C.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 6.2082 6.2082 1 29 14.1685 0.0007561 ***
## Site 2.7233 1.3616 2 12 3.1076 0.0817498 .
## Condition 1.6409 1.6409 1 14 3.7450 0.0734301 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(MC.symb.d13C.mod.FINAL)
posthoc<-emmeans(MC.symb.d13C.mod.FINAL, pairwise~Period)
CLD(posthoc, Letters=letters)## Period emmean SE df lower.CL upper.CL .group
## Bleaching -15.9 0.223 16.3 -16.3 -15.4 a
## Recovery -15.2 0.223 16.3 -15.7 -14.7 b
##
## Results are averaged over the levels of: Site, Condition
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(MC.symb.d13C.mod.FINAL), ylab="Mcap symb d13C")
M. capitata, Carbon isotope, host - symbiont: δ13CH-S
######################
# d13C.H.S.
######################
# notes on what difference in H-S isotope means
# *if H-S is (+)....*
# ex: host = -15 d^13^C and symbiont = -16 d^13^C ... difference = +1
# - the host is enriched in ^13^C, less negative and "heavier"
# *if H-S is (-)....*
# ex: host = -16 d^13^C and symbiont = -15 d^13^C ... difference = -1
# the host is depleted in ^13^C, more negative and "lighter"
# *with greater heterotrophy...*
# expect larger difference in H:S, expect host to be more depleted (more negative relative to symbiont)
# expect more (-) H/S
full.test<-lmer(d13C.H.S~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
anova(full, type=2) # Period and Site effect
#### compare models ####
# full model
full<-lmer(d13C.H.S~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(d13C.H.S~Period+Site+Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
MC.d13C.H.S.mod.FINAL<-lmer(d13C.H.S~Period+Site+Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(MC.d13C.H.S.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 2.75204 2.75204 1 29.010 12.7053 0.001286 **
## Site 0.45246 0.22623 2 12.001 1.0444 0.381797
## Condition 0.99585 0.99585 1 14.006 4.5975 0.050048 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(MC.d13C.H.S.mod.FINAL)
posthoc<-emmeans(MC.d13C.H.S.mod.FINAL, pairwise~Period)
CLD(posthoc, Letters=letters)## Period emmean SE df lower.CL upper.CL .group
## Recovery -0.112 0.102 25.1 -0.322 0.0991 a
## Bleaching 0.317 0.102 25.1 0.106 0.5274 b
##
## Results are averaged over the levels of: Site, Condition
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05posthoc<-emmeans(MC.d13C.H.S.mod.FINAL, pairwise~Condition)
CLD(posthoc, Letters=letters)## Condition emmean SE df lower.CL upper.CL .group
## NB -0.030 0.103 23 -0.2438 0.184 a
## B 0.235 0.103 23 0.0212 0.449 a
##
## Results are averaged over the levels of: Period, Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(MC.d13C.H.S.mod.FINAL), ylab="Mcap d13C H-S")
M. capitata, Nitrogen isotope, host: δ15NH
#############
# host d15N
#############
full.test<-lmer(host.d15N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
anova(full, type=2) # only Site effect
#### compare models ####
# full model
full<-lmer(host.d15N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(host.d15N~Period+Site+Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
MC.host.d15N.mod.FINAL<-lmer(host.d15N~Period+Site+Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(MC.host.d15N.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 0.0022230 0.0022230 1 29.003 0.7527 0.39275
## Site 0.0245142 0.0122571 2 12.000 4.1500 0.04267 *
## Condition 0.0000204 0.0000204 1 14.001 0.0069 0.93495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(MC.host.d15N.mod.FINAL)
posthoc<-emmeans(MC.host.d15N.mod.FINAL, pairwise~Site)
CLD(posthoc, Letters=letters)## Site emmean SE df lower.CL upper.CL .group
## Reef 25 1.42 0.0173 12 1.38 1.46 a
## Reef 44 1.45 0.0173 12 1.41 1.48 ab
## HIMB 1.49 0.0173 12 1.45 1.53 b
##
## Results are averaged over the levels of: Period, Condition
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05plot(allEffects(MC.host.d15N.mod.FINAL), ylab="Mcap host d15N")
M. capitata, Nitrogen isotope, symbiont: δ15NS
#############
# symb d15N
#############
full.test<-lmer(symb.d15N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
anova(full, type=2) # only Site effect
#### compare models ####
# full model
full<-lmer(symb.d15N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(symb.d15N~Period+Site+Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
MC.symb.d15N.mod.FINAL<-lmer(symb.d15N~Period+Site+Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(MC.symb.d15N.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 0.0045596 0.0045596 1 43.012 0.7846 0.3807
## Site 0.0265316 0.0132658 2 12.011 2.2828 0.1445
## Condition 0.0022034 0.0022034 1 43.012 0.3792 0.5413# rand(MC.symb.d15N.mod.FINAL)
plot(allEffects(MC.symb.d15N.mod.FINAL), ylab="Mcap symb d15N")
M. capitata, Nitrogen isotope, host - symbiont: δ15NH-S
##################################
############ host-symb ###########
######################
# d15N host-symb
######################
full.test<-lmer(d15N.H.S~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
anova(full, type=2) # Period and Site effect
#### compare models ####
# full model
full<-lmer(d15N.H.S~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(d15N.H.S~Period+Site+Condition+Site:Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
MC.d15N.H.S.mod.FINAL<-lmer(d15N.H.S~Period+Site+Condition+Site:Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(MC.d15N.H.S.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 1.6007 1.60067 1 29 1.7745 0.19320
## Site 0.8085 0.40423 2 24 0.4481 0.64406
## Condition 0.3954 0.39544 1 24 0.4384 0.51421
## Site:Condition 5.4921 2.74604 2 24 3.0442 0.06633 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(MC.d15N.H.S.mod.FINAL)
plot(allEffects(MC.d15N.H.S.mod.FINAL), ylab="Mcap d15N H-S")
M. capitata, C:N, host: C:NH
#############
# host C:N
#############
full.test<-lmer(host.C.N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
anova(full, type=2) # only Site effect
#### compare models ####
# full model
full<-lmer(host.C.N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(host.C.N~Period+Site+Condition+Period:Site+Period:Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
MC.host.CN.mod.FINAL<-lmer(host.C.N~Period+Site+Condition+Period:Site+Period:Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(MC.host.CN.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 3.1944 3.1944 1 26 24.6099 3.726e-05 ***
## Site 0.3750 0.1875 2 12 1.4445 0.27410
## Condition 0.4707 0.4707 1 14 3.6265 0.07762 .
## Period:Site 0.8536 0.4268 2 26 3.2881 0.05333 .
## Period:Condition 0.5681 0.5681 1 26 4.3769 0.04634 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(MC.host.CN.mod.FINAL)
posthoc<-emmeans(MC.host.CN.mod.FINAL, pairwise~Period)
CLD(posthoc, Letters=letters)## Period emmean SE df lower.CL upper.CL .group
## Bleaching 7.05 0.0885 21.5 6.86 7.23 a
## Recovery 7.51 0.0885 21.5 7.32 7.69 b
##
## Results are averaged over the levels of: Site, Condition
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05posthoc<-emmeans(MC.host.CN.mod.FINAL, pairwise~Condition|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Condition emmean SE df lower.CL upper.CL .group
## NB 7.04 0.115 42.1 6.80 7.27 a
## B 7.06 0.115 42.1 6.82 7.29 a
##
## Period = Recovery:
## Condition emmean SE df lower.CL upper.CL .group
## NB 7.30 0.115 42.1 7.07 7.53 a
## B 7.71 0.115 42.1 7.48 7.94 b
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(MC.host.CN.mod.FINAL), ylab="Mcap host C:N")
M. capitata, C:N, symb: C:NS
#############
# symb C:N
#############
full.test<-lmer(symb.C.N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
anova(full, type=2) # Period and Site effect
#### compare models ####
# full model
full<-lmer(symb.C.N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(symb.C.N~Period+Site+Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=MC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
MC.symb.C.N.mod.FINAL<-lmer(symb.C.N~Period+Site+Condition+(1|Sample.ID)+(1|Pair), data=MC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(MC.symb.C.N.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 0.053403 0.053403 1 42.552 3.7251 0.06028 .
## Site 0.047493 0.023746 2 11.989 1.6564 0.23164
## Condition 0.014074 0.014074 1 41.094 0.9817 0.32757
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(MC.symb.C.N.mod.FINAL)
plot(allEffects(MC.symb.C.N.mod.FINAL), ylab="Mcap symb C:N")
Porites compressa
Linear Mixed Effect Models
Results from these models are summarized in Table 1 and in Table S7
# PC data only
PC.iso.data<-iso.df[(iso.df$Species=="PC"),]
PC.iso.tests<-PC.iso.data # for transformations
## tests for isotopes all data
for(i in c(10:17)){
Y<-PC.iso.tests[,i]
full<-lmer(Y~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
R <- resid(full) #save glm residuals
Y.shapiro <- shapiro.test(R) #runs a normality test on residuals
print(Y.shapiro) # null = normally distrubuted (P<0.05 = non-normal
op<-par(mfrow = c(2,2), mar=c(5,4,1,2), pty="sq")
plot(full, add.smooth = FALSE, which=1)
QQ <- qqnorm(R, main = colnames(PC.iso.tests)[i])
QQline <- qqline(R)
hist(R, xlab="Residuals", main = colnames(PC.iso.tests)[i])
}P. compressa, Carbon isotope, host: δ13CH
#############
# host d13C
#############
full.test<-lmer(host.d13C~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
anova(full.test, type=2) # 3 way interaction, use full model# final model
PC.host.d13C.mod.FINAL<-lmer(host.d13C~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(PC.host.d13C.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 0.02817 0.02817 1 24 0.1110 0.741858
## Site 0.13633 0.06816 2 12 0.2687 0.768843
## Condition 1.29143 1.29143 1 12 5.0911 0.043497 *
## Period:Site 0.32033 0.16017 2 24 0.6314 0.540447
## Period:Condition 2.20417 2.20417 1 24 8.6893 0.007023 **
## Site:Condition 1.44558 0.72279 2 12 2.8494 0.097148 .
## Period:Site:Condition 2.01433 1.00717 2 24 3.9705 0.032387 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(PC.host.d13C.mod.FINAL)
posthoc<-emmeans(PC.host.d13C.mod.FINAL, pairwise~Condition|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Condition emmean SE df lower.CL upper.CL .group
## NB -15.6 0.239 24.7 -16.1 -15.1 a
## B -15.5 0.239 24.7 -16.0 -15.0 a
##
## Period = Recovery:
## Condition emmean SE df lower.CL upper.CL .group
## NB -15.9 0.239 24.7 -16.4 -15.4 a
## B -15.1 0.239 24.7 -15.5 -14.6 b
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05posthoc<-emmeans(PC.host.d13C.mod.FINAL, pairwise~Condition|Period/Site)
CLD(posthoc, Letters=letters)## Period = Bleaching, Site = HIMB:
## Condition emmean SE df lower.CL upper.CL .group
## NB -15.8 0.415 24.7 -16.7 -14.9 a
## B -15.6 0.415 24.7 -16.4 -14.7 a
##
## Period = Bleaching, Site = Reef 25:
## Condition emmean SE df lower.CL upper.CL .group
## B -15.3 0.415 24.7 -16.2 -14.5 a
## NB -15.2 0.415 24.7 -16.0 -14.3 a
##
## Period = Bleaching, Site = Reef 44:
## Condition emmean SE df lower.CL upper.CL .group
## NB -15.7 0.415 24.7 -16.6 -14.9 a
## B -15.6 0.415 24.7 -16.4 -14.7 a
##
## Period = Recovery, Site = HIMB:
## Condition emmean SE df lower.CL upper.CL .group
## NB -16.5 0.415 24.7 -17.3 -15.6 a
## B -14.4 0.415 24.7 -15.3 -13.6 b
##
## Period = Recovery, Site = Reef 25:
## Condition emmean SE df lower.CL upper.CL .group
## NB -15.4 0.415 24.7 -16.2 -14.5 a
## B -15.3 0.415 24.7 -16.2 -14.5 a
##
## Period = Recovery, Site = Reef 44:
## Condition emmean SE df lower.CL upper.CL .group
## NB -15.9 0.415 24.7 -16.7 -15.0 a
## B -15.4 0.415 24.7 -16.3 -14.5 a
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(PC.host.d13C.mod.FINAL), ylab="Pcom host.d13C")
P. compressa, Carbon isotope, symbiont: δ13CS
#############
# symb d13C
#############
full.test<-lmer(symb.d13C~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
anova(full.test, type=2) # no effects
#### compare models ####
# full model
full<-lmer(symb.d13C~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(symb.d13C~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
PC.symb.d13C.mod.FINAL<-lmer(symb.d13C~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(PC.symb.d13C.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 1.0140 1.01400 1 28 1.6317 0.21196
## Site 0.1797 0.08985 2 12 0.1446 0.86686
## Condition 2.2939 2.29386 1 14 3.6913 0.07530 .
## Period:Condition 2.6460 2.64600 1 28 4.2579 0.04845 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(PC.symb.d13C.mod.FINAL)
posthoc<-emmeans(PC.symb.d13C.mod.FINAL, pairwise~Condition|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Condition emmean SE df lower.CL upper.CL .group
## NB -15.7 0.28 39.3 -16.2 -15.1 a
## B -15.5 0.28 39.3 -16.1 -15.0 a
##
## Period = Recovery:
## Condition emmean SE df lower.CL upper.CL .group
## NB -15.8 0.28 39.3 -16.4 -15.3 a
## B -14.9 0.28 39.3 -15.4 -14.3 b
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(PC.symb.d13C.mod.FINAL), ylab="Pcom symb d13C")
P. compressa, Carbon isotope, host - symbiont: δ13CH-S
######################
# d13C.H.S.
###################### no effects
full.test<-lmer(d13C.H.S~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
anova(full.test, type=2)
#### compare models ####
# full model
full<-lmer(d13C.H.S~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(d13C.H.S~Period+Site+Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
PC.d13C.H.S.mod.FINAL<-lmer(d13C.H.S~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(PC.d13C.H.S.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 0.69338 0.69338 1 54 2.2867 0.1363
## Site 0.52908 0.26454 2 54 0.8724 0.4237
## Condition 0.12604 0.12604 1 54 0.4157 0.5218
## Period:Condition 0.01504 0.01504 1 54 0.0496 0.8246# rand(PC.d13C.H.S.mod.FINAL)
plot(allEffects(PC.d13C.H.S.mod.FINAL), ylab="Pcom symb d13C.H-S")
P. compressa, Nitrogen isotope, host: δ15NH
#############
# host d15N
############# Period, Site, Period:Condition effects
full.test<-lmer(host.d15N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
anova(full.test, type=2)
#### compare models ####
# full model
full<-lmer(host.d15N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(host.d15N~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
PC.host.d15N.mod.FINAL<-lmer(host.d15N~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(PC.host.d15N.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 2.7735 2.7735 1 28.000 6.7953 0.014485 *
## Site 8.4827 4.2413 2 11.961 10.3917 0.002423 **
## Condition 0.0985 0.0985 1 13.955 0.2413 0.630923
## Period:Condition 2.0535 2.0535 1 28.000 5.0313 0.032991 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(PC.host.d15N.mod.FINAL)
posthoc<-emmeans(PC.host.d15N.mod.FINAL, pairwise~Period)
CLD(posthoc, Letters=letters)## Period emmean SE df lower.CL upper.CL .group
## Bleaching 4.91 0.129 28.4 4.65 5.17 a
## Recovery 5.34 0.129 28.4 5.08 5.60 b
##
## Results are averaged over the levels of: Site, Condition
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05posthoc<-emmeans(PC.host.d15N.mod.FINAL, pairwise~Site)
CLD(posthoc, Letters=letters)## Site emmean SE df lower.CL upper.CL .group
## Reef 25 4.76 0.172 12 4.39 5.14 a
## Reef 44 4.84 0.172 12 4.47 5.22 a
## HIMB 5.76 0.172 12 5.39 6.14 b
##
## Results are averaged over the levels of: Period, Condition
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05posthoc<-emmeans(PC.host.d15N.mod.FINAL, pairwise~Condition|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Condition emmean SE df lower.CL upper.CL .group
## B 4.77 0.182 51.3 4.41 5.14 a
## NB 5.05 0.182 51.3 4.68 5.41 a
##
## Period = Recovery:
## Condition emmean SE df lower.CL upper.CL .group
## NB 5.11 0.182 51.3 4.74 5.47 a
## B 5.57 0.182 51.3 5.21 5.94 a
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(PC.host.d15N.mod.FINAL), ylab="Pcom host.d15N")
P. compressa, Nitrogen isotope, symbiont: δ15NS
#############
# symb d15N
############# Site, Period:Condition effects
full.test<-lmer(symb.d15N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
anova(full.test, type=2) # no effects
#### compare models ####
# full model
full<-lmer(symb.d15N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(symb.d15N~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
PC.symb.d15N.mod.FINAL<-lmer(symb.d15N~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(PC.symb.d15N.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 2.4402 2.4402 1 28.000 2.6549 0.114435
## Site 9.4572 4.7286 2 12.011 5.1446 0.024329 *
## Condition 1.4517 1.4517 1 14.016 1.5794 0.229385
## Period:Condition 7.1415 7.1415 1 28.000 7.7698 0.009436 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(PC.symb.d15N.mod.FINAL)
posthoc<-emmeans(PC.symb.d15N.mod.FINAL, pairwise~Site)
CLD(posthoc, Letters=letters)## Site emmean SE df lower.CL upper.CL .group
## Reef 44 3.75 0.272 12 3.16 4.35 a
## Reef 25 3.98 0.272 12 3.38 4.57 ab
## HIMB 4.92 0.272 12 4.32 5.51 b
##
## Results are averaged over the levels of: Period, Condition
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05posthoc<-emmeans(PC.symb.d15N.mod.FINAL, pairwise~Condition|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Condition emmean SE df lower.CL upper.CL .group
## NB 3.88 0.28 50.3 3.32 4.44 a
## B 4.95 0.28 50.3 4.39 5.52 b
##
## Period = Recovery:
## Condition emmean SE df lower.CL upper.CL .group
## B 3.86 0.28 50.3 3.30 4.42 a
## NB 4.17 0.28 50.3 3.60 4.73 a
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(PC.symb.d15N.mod.FINAL), ylab="Pcom symb d15N")
P. compressa, Nitrogen isotope, host - symbiont: δ15NH-S
######################
# d15N host-symb
######################
full.test<-lmer(d15N.H.S~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
anova(full.test, type=2)
#### compare models ####
# full model
full<-lmer(d15N.H.S~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(d15N.H.S~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
PC.d15N.H.S.mod.FINAL<-lmer(d15N.H.S~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(PC.d15N.H.S.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 10.2920 10.2920 1 27.996 7.8431 0.009144 **
## Site 0.6065 0.3032 2 11.998 0.2311 0.797124
## Condition 0.8763 0.8763 1 13.996 0.6678 0.427511
## Period:Condition 16.4850 16.4850 1 27.996 12.5625 0.001405 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(PC.d15N.H.S.mod.FINAL)
posthoc<-emmeans(PC.d15N.H.S.mod.FINAL, pairwise~Condition|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Condition emmean SE df lower.CL upper.CL .group
## B -0.173 0.335 49.9 -0.846 0.499 a
## NB 1.167 0.335 49.9 0.494 1.839 b
##
## Period = Recovery:
## Condition emmean SE df lower.CL upper.CL .group
## NB 0.947 0.335 49.9 0.274 1.619 a
## B 1.703 0.335 49.9 1.031 2.376 a
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(PC.d15N.H.S.mod.FINAL), ylab="Pcom symb d15N.H-S")
P. compressa, C:N, host: C:NH
#############
# host C:N
#############
full.test<-lmer(host.C.N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
anova(full.test, type=2) # no effects
#### compare models ####
# full model
full<-lmer(host.C.N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(host.C.N~Period+Site+Condition+Period:Site+Period:Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
PC.host.C.N.mod.FINAL<-lmer(host.C.N~Period+Site+Condition+Period:Site+Period:Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(PC.host.C.N.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 5.8737 5.8737 1 26.000 33.8692 3.924e-06 ***
## Site 0.1503 0.0752 2 11.999 0.4334 0.6580649
## Condition 0.0396 0.0396 1 13.999 0.2284 0.6400871
## Period:Site 2.3244 1.1622 2 26.000 6.7016 0.0044953 **
## Period:Condition 3.4696 3.4696 1 26.000 20.0067 0.0001351 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rand(PC.host.C.N.mod.FINAL)
posthoc<-emmeans(PC.host.C.N.mod.FINAL, pairwise~Period)
CLD(posthoc, Letters=letters)## Period emmean SE df lower.CL upper.CL .group
## Bleaching 6.29 0.0963 23.1 6.09 6.49 a
## Recovery 6.92 0.0963 23.1 6.72 7.11 b
##
## Results are averaged over the levels of: Site, Condition
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05posthoc<-emmeans(PC.host.C.N.mod.FINAL, pairwise~Site|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Site emmean SE df lower.CL upper.CL .group
## HIMB 6.15 0.167 23.1 5.81 6.50 a
## Reef 25 6.26 0.167 23.1 5.91 6.60 a
## Reef 44 6.46 0.167 23.1 6.12 6.81 a
##
## Period = Recovery:
## Site emmean SE df lower.CL upper.CL .group
## Reef 44 6.54 0.167 23.1 6.19 6.88 a
## Reef 25 7.09 0.167 23.1 6.75 7.44 a
## HIMB 7.12 0.167 23.1 6.77 7.46 a
##
## Results are averaged over the levels of: Condition
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05posthoc<-emmeans(PC.host.C.N.mod.FINAL, pairwise~Condition|Period)
CLD(posthoc, Letters=letters)## Period = Bleaching:
## Condition emmean SE df lower.CL upper.CL .group
## B 6.09 0.134 45.5 5.82 6.36 a
## NB 6.49 0.134 45.5 6.22 6.76 b
##
## Period = Recovery:
## Condition emmean SE df lower.CL upper.CL .group
## NB 6.64 0.134 45.5 6.37 6.91 a
## B 7.19 0.134 45.5 6.92 7.46 b
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(PC.host.C.N.mod.FINAL), ylab="Pcom host.C.N")
P. compressa, C:N, symbiont: C:NS
#############
# symb CN
#############
full.test<-lmer(symb.C.N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
anova(full.test, type=2) # no effects
#### compare models ####
# full model
full<-lmer(symb.C.N~Period*Site*Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
# additive model
add<-lmer(symb.C.N~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), REML=FALSE, data=PC.iso.tests, na.action=na.exclude)
anova(full, add) # no difference# final model
PC.symb.C.N.mod.FINAL<-lmer(symb.C.N~Period+Site+Condition+Period:Condition+(1|Sample.ID)+(1|Pair), data=PC.iso.tests, na.action=na.exclude)
# ANOVA table
anova(PC.symb.C.N.mod.FINAL, type=2)## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 3.2025 3.2025 1 5.5729 2.5883 0.1625
## Site 0.2094 0.1047 2 2.0000 0.0846 0.9220
## Condition 3.8472 3.8472 1 5.5729 3.1094 0.1321
## Period:Condition 2.0104 2.0104 1 5.5729 1.6248 0.2530# rand(PC.symb.C.N.mod.FINAL)
plot(allEffects(PC.symb.C.N.mod.FINAL), ylab="Pcom symb C:N")
Isotope figures
Figures here are for: δ13CHost, δ13CSymb, δ15NHost, δ15NSymb, δ13CHost-Symb, δ15NHost-Symb, C:NHost, C:NSymb.
These figures together make up Figure 7 (isotope data, both species). Isotope figures are mean +/-SE with Montipora capitata on the left, Porites compressa on the right.
# isotope summaries by species
MC.iso<-iso.summary[(iso.summary$Species=='MC'),] # MC isotopes
PC.iso<-iso.summary[(iso.summary$Species=='PC'),] # PC isotopes
MC.iso$Site<-factor(MC.iso$Site, levels=c("Reef 44", "Reef 25", "HIMB"))
# general formatting for all isotope figures
Fig.formatting<-(theme_classic())+
theme(text=element_text(size=8),
axis.line=element_blank(),
legend.justification=c(1,1), legend.position=c(1,1),
legend.background = element_rect("transparent", colour=NA),
legend.text.align = 0,
legend.text=element_text(size=8),
legend.title = element_blank(),
panel.border = element_rect(fill=NA, colour = "black", size=1),
aspect.ratio=.7,
axis.ticks.length=unit(0.25, "cm"),
axis.text.y=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=8),
axis.text.x=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=5)) +
theme(legend.key.size = unit(0.4, "cm")) +
theme(panel.spacing=unit(c(0, 0, 0, 0), "cm"))
pd <- position_dodge(0.71) #offset for error bars and columns
################### d13C Host
Fig.MC.iso.host.d13C<-ggplot(MC.iso, aes(x=Period:Site, y=host.d13C, fill=Condition, group=Condition)) +
geom_errorbar(aes(ymin=host.d13C-host.d13C.SE, ymax=host.d13C+host.d13C.SE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))) +
geom_point(aes(fill=Condition), position=pd, size=3, pch=21) +
scale_y_continuous(limits=c(-19, -12)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) +
geom_vline(xintercept = 3.5, linetype="dotted")+
Fig.formatting
################### d13C Symbiont
Fig.MC.iso.symb.d13C<-ggplot(MC.iso, aes(x=Period:Site, y=symb.d13C, fill=Condition, group=Condition)) +
geom_errorbar(aes(ymin=symb.d13C-symb.d13C.SE, ymax=symb.d13C+symb.d13C.SE),size=.5, width=0, position=pd) +
xlab("") +
geom_point(aes(fill=Condition), position=pd, size=3, pch=21) +
scale_y_continuous(limits=c(-19, -12)) +
ylab(expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)"))) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) + guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted")+ Fig.formatting
################### d13C Host-Symbiont
Fig.MC.iso.d13C.H.S<-ggplot(MC.iso, aes(x=Period:Site, y=d13C.H.S, fill=Condition, group=Condition)) +
geom_point(stat="identity", position = pd) +
geom_hline(yintercept=0, linetype="solid", color = "red") +
geom_errorbar(aes(ymin=d13C.H.S-d13C.H.S.SE, ymax=d13C.H.S+d13C.H.S.SE),size=.5, width=0, position=pd) +
xlab("") +
geom_point(aes(fill=Condition), position=pd, size=3, pch=21) +
scale_y_continuous(limits=c(-1.5, 1.5)) +
ylab(expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)"))) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) + guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted") +
Fig.formatting
################### d15N Host
Fig.MC.iso.host.d15N<-ggplot(MC.iso, aes(x=Period:Site, y=host.d15N, fill=Condition, group=Condition)) +
geom_errorbar(aes(ymin=host.d15N-host.d15N.SE, ymax=host.d15N+host.d15N.SE),size=.5, width=0, position=pd) +
xlab("") +
geom_point(aes(fill=Condition), position=pd, size=3, pch=21) +
scale_y_continuous(limits=c(2, 8)) +
ylab(expression(paste(delta^{15}, N[H], " (\u2030, air)"))) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached)))+ guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted") +
Fig.formatting
################### d15N Symbiont
Fig.MC.iso.symb.d15N<-ggplot(MC.iso, aes(x=Period:Site, y=symb.d15N, fill=Condition, group=Condition)) +
geom_errorbar(aes(ymin=symb.d15N-symb.d15N.SE, ymax=symb.d15N+symb.d15N.SE),size=.5, width=0, position=pd) +
xlab("") +
geom_point(aes(fill=Condition), position=pd, size=3, pch=21) +
scale_y_continuous(limits=c(2, 8)) +
ylab(expression(paste(delta^{15}, N[S], " (\u2030, air)"))) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached)))+ guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted")+
Fig.formatting
################### d15N Host-Symbiont
Fig.MC.iso.d15N.H.S<-ggplot(MC.iso, aes(x=Period:Site, y=d15N.H.S, fill=Condition, group=Condition)) +
geom_point(stat="identity", position = pd) +
geom_hline(yintercept=0, linetype="solid", color = "red") +
geom_errorbar(aes(ymin=d15N.H.S-d15N.H.S.SE, ymax=d15N.H.S+d15N.H.S.SE),size=.5, width=0, position=pd) +
xlab("") +
geom_point(aes(fill=Condition), position=pd, size=3, pch=21) +
scale_y_continuous(limits=c(-3, 5)) +
ylab(expression(paste(delta^{15}, N[H-S], " (\u2030, air)"))) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) + guides(fill="none")+
geom_vline(xintercept = 3.5, linetype="dotted")+
Fig.formatting
################### C:N Host
Fig.MC.iso.hostCN<-ggplot(MC.iso, aes(x=Period:Site, y=host.C.N, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=host.C.N-host.C.N.SE, ymax=host.C.N+host.C.N.SE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("C:N "[H]))) +
coord_cartesian(ylim=c(5,10)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached)))+
geom_vline(xintercept = 3.5, linetype="dotted")+
Fig.formatting
################### C:N Symbiont
Fig.MC.iso.symbCN<-ggplot(MC.iso, aes(x=Period:Site, y=symb.C.N, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=symb.C.N-symb.C.N.SE, ymax=symb.C.N+symb.C.N.SE),size=.5, width=0, position=pd) +
xlab("") +
coord_cartesian(ylim=c(5,10)) +
scale_x_discrete(labels=Period.Reef) +
ylab(expression(paste("C:N "[S]))) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) + guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted") +
Fig.formatting
#####
###################
#############################
###################
##### Porites isotopes
pd<-position_dodge(0.71) #offset for error bars
PC.iso$Site<-factor(PC.iso$Site, levels=c("Reef 44", "Reef 25", "HIMB"))
################### d13C Host
Fig.PC.iso.host.d13C<-ggplot(PC.iso, aes(x=Period:Site, y=host.d13C, fill=Condition, group=Condition)) +
geom_errorbar(aes(ymin=host.d13C-host.d13C.SE, ymax=host.d13C+host.d13C.SE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))) +
geom_point(aes(fill=Condition), position=pd, size=3, pch=21) +
scale_y_continuous(limits=c(-19, -12)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached)))+ guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted")+
Fig.formatting
################### d13C Symbiont
Fig.PC.iso.symb.d13C<-ggplot(PC.iso, aes(x=Period:Site, y=symb.d13C, fill=Condition, group=Condition)) +
geom_errorbar(aes(ymin=symb.d13C-symb.d13C.SE, ymax=symb.d13C+symb.d13C.SE),size=.5, width=0, position=pd) +
xlab("") +
geom_point(aes(fill=Condition), position=pd, size=3, pch=21) +
scale_y_continuous(limits=c(-19, -12)) +
ylab(expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)"))) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) + guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted")+ Fig.formatting
################### d13C Host-Symbiont
Fig.PC.iso.d13C.H.S<-ggplot(PC.iso, aes(x=Period:Site, y=d13C.H.S, fill=Condition, group=Condition)) +
geom_point(stat="identity", position = pd) +
geom_hline(yintercept=0, linetype="solid", color = "red") +
geom_errorbar(aes(ymin=d13C.H.S-d13C.H.S.SE, ymax=d13C.H.S+d13C.H.S.SE),size=.5, width=0, position=pd) +
xlab("") +
geom_point(aes(fill=Condition), position=pd, size=3, pch=21) +
scale_y_continuous(limits=c(-1.5, 1.5)) +
ylab(expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)"))) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) + guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted") +
Fig.formatting
################### d15N Host
Fig.PC.iso.host.d15N<-ggplot(PC.iso, aes(x=Period:Site, y=host.d15N, fill=Condition, group=Condition)) +
geom_errorbar(aes(ymin=host.d15N-host.d15N.SE, ymax=host.d15N+host.d15N.SE),size=.5, width=0, position=pd) +
xlab("") +
geom_point(aes(fill=Condition), position=pd, size=3, pch=21) +
scale_y_continuous(limits=c(2, 8)) +
ylab(expression(paste(delta^{15}, N[H], " (\u2030, air)"))) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) + guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted") +
Fig.formatting
################### d15N Symbiont
Fig.PC.iso.symb.d15N<-ggplot(PC.iso, aes(x=Period:Site, y=symb.d15N, fill=Condition, group=Condition)) +
geom_errorbar(aes(ymin=symb.d15N-symb.d15N.SE, ymax=symb.d15N+symb.d15N.SE),size=.5, width=0, position=pd) +
xlab("") +
geom_point(aes(fill=Condition), position=pd, size=3, pch=21) +
scale_y_continuous(limits=c(2, 8)) +
ylab(expression(paste(delta^{15}, N[S], " (\u2030, air)"))) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached)))+ guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted")+
Fig.formatting
################### d15N Host-Symbiont
Fig.PC.iso.d15N.H.S<-ggplot(PC.iso, aes(x=Period:Site, y=d15N.H.S, fill=Condition, group=Condition)) +
geom_point(stat="identity", position = pd) +
geom_hline(yintercept=0, linetype="solid", color = "red") +
geom_errorbar(aes(ymin=d15N.H.S-d15N.H.S.SE, ymax=d15N.H.S+d15N.H.S.SE),size=.5, width=0, position=pd) +
xlab("") +
geom_point(aes(fill=Condition), position=pd, size=3, pch=21) +
scale_y_continuous(limits=c(-3, 5)) +
ylab(expression(paste(delta^{15}, N[H-S], " (\u2030, air)"))) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) + guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted") +
Fig.formatting
################### C:N Host
Fig.PC.iso.hostCN<-ggplot(PC.iso, aes(x=Period:Site, y=host.C.N, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=host.C.N-host.C.N.SE, ymax=host.C.N+host.C.N.SE),size=.5, width=0, position=pd) +
xlab("") +
ylab(expression(paste("C:N "[H]))) +
coord_cartesian(ylim=c(5,10)) +
scale_x_discrete(labels=Period.Reef) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) + guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted")+
Fig.formatting
################### C:N Symbiont
Fig.PC.iso.symbCN<-ggplot(PC.iso, aes(x=Period:Site, y=symb.C.N, fill=Condition, group=Condition)) +
geom_bar(colour="black", stat="identity", position = pd, width=0.7) +
geom_errorbar(aes(ymin=symb.C.N-symb.C.N.SE, ymax=symb.C.N+symb.C.N.SE),size=.5, width=0, position=pd) +
xlab("") +
coord_cartesian(ylim=c(5,10)) +
scale_x_discrete(labels=Period.Reef) +
ylab(expression(paste("C:N "[S]))) +
scale_fill_manual(values=c('gray92','black'),
labels=rep(expression(Bleached, Non-Bleached))) + guides(fill="none") +
geom_vline(xintercept = 3.5, linetype="dotted") +
Fig.formatting
################
################ combine figures, as in the manuscript
################
plot_grid(Fig.MC.iso.host.d13C, Fig.PC.iso.host.d13C,
Fig.MC.iso.symb.d13C, Fig.PC.iso.symb.d13C,
labels = c("a", "b", "c", "d"), ncol = 2, axis="lr")
plot_grid(Fig.MC.iso.d13C.H.S, Fig.PC.iso.d13C.H.S,
Fig.MC.iso.host.d15N, Fig.PC.iso.host.d15N,
labels = c("e", "f", "g", "h"), ncol = 2, axis="lr")
plot_grid(Fig.MC.iso.symb.d15N, Fig.PC.iso.symb.d15N,
Fig.MC.iso.d15N.H.S, Fig.PC.iso.d15N.H.S,
labels = c("i", "j", "k", "l"), ncol = 2, axis="lr")
Next batch of figures are for molar concentrations of carbon and nitrogen in host and symbionts.
Figure S2 (C:N host and symbiont, both species). C:N molar ratios are mean +/-SE with Montipora capitata on the left, Porites compressa on the right.
plot_grid(Fig.MC.iso.hostCN, Fig.PC.iso.hostCN,
Fig.MC.iso.symbCN, Fig.PC.iso.symbCN,
labels = c("a", "b", "c", "d"), ncol = 2, axis="lr")
Isotope mass balance
An isotope mass balance was modeled to measure changes in tissue biomass composition on holobiont (host + symbiont) δ13C values during bleaching recovery, following Hayes (2001).
First, the isotopic composition of the holobiont (δ13CHolobiont) was modeled for each time period using ash-free dry weight (AFDW) and assuming holobiont AFDW is 95% host vs. 5% symbiont mass.
Second, the measured proportion of biomass compounds (i.e., % of proteins, lipids, carbohydrates) and δ13CHolobiont were used to estimate compound-specific isotopic values (δ13CCompound) for each compound. The influence of changes in biomass composition on δ13CHolobiont values (i.e., observed-δ13CHolobiont) during bleaching recovery δ13CCompound estimates were applied to the same colonies in January 2015 to determine expected-δ13CHolobiont. The relationship between observed and expected δ13CHolobiont was evaluated using a linear regression.
Figures displayed here can be found in Figure S3 (δ13CCompound) and Figure 8 (observed- vs. expected-δ13CHolobiont)
######## Mass balance attempt#######
# δBiomass = XCNA*δNA + XCProt*δProt +XCSacc*δSacc + XCLip*δLip
# XC should sum to 1 if all masses are measured correctly, refers to mole fractions of carbon
# δ terms refer to mass-weighted average isotopic compositions for all compounds within the indicated classe
Mass.df<-data.full[, !names(data.full) %in% c("symb.comm", "chla", "chlc2", "chltot", "energy.gdw", "host.ug.N", "host.d15N", "host.ug.C", "host.C.N", "symb.ug.N", "symb.d15N", "symb.ug.C", "symb.C.N", "d15N.H.S", "d13C.H.S")]
# mass proportions for formulas below
h.B<-0.95 # bleached hosts
s.B<-0.05 # bleached symb
h.NB<-0.95 # non-bleaced hosts
s.NB<- 0.05 # non-bleached symb
#### calculate whole tissue d13C for symb + host based on literature % masses in whole tissues (above)
# assuming a 5% biomass from symbionts in bleached corals and 12 % in nonbleached (Thornhill et al. 2011)
# this is considered the 'holobiont d13C'
Mass.df<-Mass.df %>% mutate(
dC.holobiont = ifelse(Period =="Bleaching" & Condition=="B",
((h.B*host.d13C)+(s.B*symb.d13C)),
((h.NB*host.d13C)+(s.NB*symb.d13C))))
#### calculate fractional abundance for each tissue type relative to total amounts of tissue measured
Mass.df<-Mass.df %>% mutate(
f.lipid = lipid.gdw/(lipid.gdw+prot.gdw+carb.gdw),
f.prot = prot.gdw/(lipid.gdw+prot.gdw+carb.gdw),
f.carb = carb.gdw/(lipid.gdw+prot.gdw+carb.gdw))
# in calculations for mass balance...
# lipids are generally 3-5 in corals (at least in the Red Sea; Alamaru et al. 2009).
#### host mass delta
# From Hayes 2000
# d.carbs = d.prot +1 ‰
# d.lipid = d.prot -5 ‰
# First, determine d13C.prot
# d13C-for compounds
Mass.df<-Mass.df %>% mutate(
d13C.prot = dC.holobiont+(5*f.lipid)-(1*f.carb)/(f.lipid+f.carb+f.prot), #d13C-protein
d13C.lip = d13C.prot-5, # d13C-lipid
d13C.carb = d13C.prot+1) # d13C-carb
# different among species? Site and Species interactions matter
anova(lm(d13C.lip~Site*Species, data=Mass.df))
anova(lm(dC.holobiont~Site*Species, data=Mass.df))
### plot these values: by species, by state, for bleaching and recovery period
Mdf.MC.bl<-Mass.df[(Mass.df$Species=="MC" & Mass.df$Condition=="B"),] # for Montipora
Mdf.MC.nb<-Mass.df[(Mass.df$Species=="MC" & Mass.df$Condition=="NB"),]
Mdf.PC.bl<-Mass.df[(Mass.df$Species=="PC" & Mass.df$Condition=="B"),] # for Porites
Mdf.PC.nb<-Mass.df[(Mass.df$Species=="PC" & Mass.df$Condition=="NB"),]
# plot of the holobiont to lipids, carbs, proteins separated by species and bleaching states
# lines match model of full data
par(mfrow=c(1,1),mar=c(5,4.5,2,1))
# Montipora
plot(d13C.lip~dC.holobiont, data=Mdf.MC.bl, type="p", pch=1, cex=0.5, col='coral', ylim=c(-22, -5), xlim=c(-19, -12), ylab=expression(paste(delta^{13}, C[Compound], " (\u2030, V-PDB)")), xlab=expression(paste(delta^{13}, C[Holobiont], " (\u2030, V-PDB)")), main="Biomass compound d13C vs. Holobiont d13C", cex.main=0.8)
legend("topleft", legend=c("lipids", "proteins", "carbs"), col=c("coral", "dodgerblue", "forestgreen"), pch=19, y.intersp=1, bty="n", cex=0.9)
par(new=T) #nonbleached lipid
points(d13C.lip~dC.holobiont, data=Mdf.MC.nb, type="p", pch=16, cex=0.5, col='coral')
ablineclip(lm(d13C.lip~dC.holobiont, data=Mass.df),x1=min(Mass.df$dC.holobiont), x2=max(Mass.df$dC.holobiont), col="coral"); par(new=T)
#bleached protein
points(d13C.prot~dC.holobiont, data=Mdf.MC.bl, pch=1, cex=0.5, col='dodgerblue'); par(new=T)
#nonbleached protein
points(d13C.prot~dC.holobiont, data=Mdf.MC.nb, pch=16, cex=0.5, col='dodgerblue')
ablineclip(lm(d13C.prot~dC.holobiont, data=Mass.df), x1=min(Mass.df$dC.holobiont), x2=max(Mass.df$dC.holobiont), col="dodgerblue"); par(new=T)
#bleached carbs
points(d13C.carb~dC.holobiont, data=Mdf.MC.bl, pch=1, cex=0.5, col='forestgreen'); par(new=T)
#nobleached carbs
points(d13C.carb~dC.holobiont, data=Mdf.MC.nb, pch=16, cex=0.5, col='forestgreen')
ablineclip(lm(d13C.carb~dC.holobiont, data=Mass.df), x1=min(Mass.df$dC.holobiont), x2=max(Mass.df$dC.holobiont), col="forestgreen"); par(new=T)
# Porites
points(d13C.lip~dC.holobiont, data=Mdf.PC.bl, pch=2, cex=0.5, col='coral'); par(new=T)
points(d13C.lip~dC.holobiont, data=Mdf.PC.nb, pch=17, cex=0.5, col='coral')
par(new=T)
points(d13C.prot~dC.holobiont, data=Mdf.PC.bl, pch=2, cex=0.5, col='dodgerblue'); par(new=T)
points(d13C.prot~dC.holobiont, data=Mdf.PC.nb, pch=17, cex=0.5, col='dodgerblue'); par(new=T)
points(d13C.carb~dC.holobiont, data=Mdf.PC.bl, pch=2, cex=0.5, col='forestgreen'); par(new=T)
points(d13C.carb~dC.holobiont, data=Mdf.PC.nb, pch=17, cex=0.5, col='forestgreen')
legend("topright", legend=c("Montipora", "Porites"), col=c("black", "black"), pch=c(16, 17), y.intersp=1, cex=0.9, bty="n")
####################
####################
#### Now to compare the d13C-compound (bleaching) with a change in proporton of compounds (recovery)
Mdf.bl<-Mass.df[(Mass.df$Time.point=="2014 Oct"),] # for Bleaching
Mdf.rc<-Mass.df[(Mass.df$Time.point=="2015 Jan"),] # for Recovery
# Mdf.bl = bleaching data
# Mdf.rc = recovery data
calc.df<- Mdf.bl[, c(1:8, 16:22)] %>% mutate(
f.lip.rec = Mdf.rc$f.lipid,
f.prot.rec = Mdf.rc$f.prot,
f.carb.rec = Mdf.rc$f.carb)
# should value of bleached lipids be -3 depleted
# calc.df$d13C.lip<- ifelse(calc.df$Condition=="B", (calc.df$d13C.lip-3), calc.df$d13C.lip)
#### Calculate a theoretical d13C-holobiont, based on d13C-biomass
calc.df<- calc.df %>% mutate(
d13C.theor = ((d13C.prot*f.prot.rec) + (d13C.carb*f.carb.rec) + (d13C.lip*f.lip.rec)))
#### Now plot relationship between d13C-holobiont-recovery and d13C-biomass with new % (recovery)
MCcalc.df<-calc.df[(calc.df$Species=="MC"), ]
PCcalc.df<-calc.df[(calc.df$Species=="PC"), ]
summary(lm(d13C.theor~dC.holobiont, data=MCcalc.df)) #R2 = 88
summary(lm(d13C.theor~dC.holobiont, data=PCcalc.df)) #R2 = 56
###### plot
par(mfrow=c(1,1),mar=c(5,4.5,2,1))
plot(d13C.theor~dC.holobiont, data=MCcalc.df, type="p", pch=16, cex=0.7, col='black', ylim=c(-20, -12), xlim=c(-20, -12), ylab=expression(paste("expected ", delta^{13}, C[Holobiont], " (\u2030, V-PDB)")), xlab=expression(paste("observed ", delta^{13}, C[Holobiont], " (\u2030, V-PDB)")), main="Observed vs. Expected d13C-Holobiont", cex.main=0.8)
par(new=T)
points(d13C.theor~dC.holobiont, data=PCcalc.df, type="p", pch=17, cex=0.7, col='gray60')
ablineclip(lm(d13C.theor~dC.holobiont, data=MCcalc.df), col="black",
x1=min(MCcalc.df$dC.holobiont), x2=max(MCcalc.df$dC.holobiont), lty=1)
ablineclip(lm(d13C.theor~dC.holobiont, data=PCcalc.df), col="gray60",
x1=min(PCcalc.df$dC.holobiont), x2=max(PCcalc.df$dC.holobiont), lty=2)
legend("topleft", legend=c("M. capitata", "P. compressa"), col=c("black", "gray60"), pch=16, lty=c(1,2), y.intersp=0.9, cex=0.9, bty="n")
text(-15, -20, "Mcap: R2=0.88, Pcom: R2=0.56, p<0.001", cex=0.8) #add in model output
############